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ABSTRACT 



Context. We present EvoL, the new release of the Padova N-body code for cosmological simulations of galaxy formation and 
evolution. In this paper, the basic Tree + SPH code is presented and analysed, together with an overview on the software 
architectures. 

Aims. EvoL is a flexible parallel Fortran95 code, specifically designed for simulations of cosmological structure formation on 
cluster, galactic and sub-galactic scales. 

Methods. EvoL is a fully Lagrangian self-adaptive code, based on the classical Oct-tree by Barnes & Hut (1986) and on the 
Smoothed Particle Hydrodynamics algorithm (SPH, Lucy 1977). It includes special features such as adaptive softening lengths 
with correcting extra-terms, and modern formulations of SPH and artificial viscosity. It is designed to be run in parallel on 
multiple CPUs to optimize the performance and save computational time. 

Results. We describe the code in detail, and present the results of a number of standard hydrodynamical tests. 
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The importance of numerical simulations in modern sci- 
ence is constantly growing, because of the complexity, the 
multi-scaling properties, and the non-linearity of many 
physical phenomena. When analytical predictions are not 
possible, we are forced to compute the evolution of physi- 
cal systems numerically. Typical examples in astrophysical 
context are the problems of structure and galaxy forma- 
tion and evolution. Over the past two decades, thanks 
to highly sophisticated cosmological and galaxy-sized nu- 
merical simulations a number of issues concerning the for- 
mation of cosmic systems and their evolution along the 
history of the Universe have been clarified. However, an 
equivalent number of new questions have been raised, so 
that complete understanding of how galaxies and clusters 
formed and evolved along the Hubble time is still out of 
our reach. This is especially true at the light of many re- 
cent observations that often appear at odds with well es- 
tablished theories (see for instance the long debated ques- 
tion of the monolithic versus hierarchical mechanism of 
galaxy formation and their many variants) and to require 



new theoretical scenarios able to match the observational 
data. 

To this aim, more and more accurate and detailed 
numerical simulations are still needed. They indeed are 
the best tool to our disposal to investigate such com- 
plex phenomena as the formation and evolution of galax- 
ies within a consistent cosmological context. A number of 
astrophysical codes for galaxy-sized and cosmological sim- 
ulations are freely and publicly available. They display a 
huge range of options and functions; each of them is best 
suited for a set of particular problems and experiments, 
and may suffer one or mo re drawbacks. Among the best 

2QQ0h . Gadget 



known, we recall Fl ash (^Frvxell et al 



fSDrinsel et alj l2QQlh and its second r elease Gadget2 
(Springel 2QQ 5j), Gasol ine (jWadslev et a l. 2004). Hydra 
(ICouchman et al.lll995|). En zo ([Norman et al.ll2007l) . and 



Vine (jWetzstein et al.ll2008f ). 
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EvoL is the new release o f the Padova N-body 
code (Pd-Tsph "Carraro^ejalJ Il998l : lUia et al.l l2QQ2l : 
Merlin fc Chiosil 12007). It is a flexible, fully Lagrangian, 
parallel, and self-adaptive N-body code, written in 
Fortran95 in a straightforward and user- friendly format. 
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EvoL describes the dynamical evolution of a system 
of interacting discrete masses (particles) moving under 
the mutual gravitational forces and/or the gravitational 
action of external bodies, plus, when appropriate, under 
mutual hydrodynamical interactions. Such particles can 
represent real discrete bodies or a fluid contained in vol- 
ume elements of suitable size. A numerical simulation of a 
physical system follows the temporal evolution of it using 
small but finite time-steps to approximate the equations of 
motion to a finite-differences problem. In the Lagrangian 
description, no reference grid is superposed to the volume 
under consideration, while the particles move under their 
mutual (or external) interactions. Each particles carries a 
mass, a position, a velocity, plus (when necessary) a list 
of physical features such as density, temperature, chemical 
composition, etc. 

To simulate the dynamical evolution of a region of 
the Universe, one has to properly model the main differ- 
ent material components, namely Dark Matter, Gas, and 
Stars), representing them with different species of parti- 
cles. Moreover, a physically robust model for the funda- 
mental interactions (in addition to gravity) is required, at 
the various scales of interest. Finally, a suitable cosmo- 
logical framework is necessary, together with a coherent 
setting of the boundary conditions and with an efficient 
algorithm to follow the temporal evolution of the system. 
EvoL is designed to respond to all of these requirements, 
improving upon the previous versions of the Padova Tree- 
SPH code under many aspects0 In the following sections, 
the main features of EvoL are presented in detail, to- 
gether with a review of some general considerations about 
the adopted algorithms whenever appropriate. 

This paper presents the basic features of the code; 
namely, the Tree-SPH (i.e. Oct-Tree plus Smoothing 
Particle Hydrodynamics) formalism, its implementation 
in the parallel architecture of the code, and the results of 
a number of classic hydrodynamic tests. The non-standard 
algorithms (e.g. radiative cooling functions, chemical evo- 
lution, star formation, energy feedback) will be presented 
in a following companion paper (Merlin et al. 2009, in 
preparation). 

The plan of the paper is as follows. In Sect. [2] we de- 
scribe the aspects of the code related to the gravitational 
interaction. In Sect. [3] we deal with the hydrodynamical 
treatment of fluids and its approximation in the SPH lan- 
guage. Sect. [4] illustrates some aspects related to the in- 
tegration technique, such as the time steps, the periodic 

^ The core of the old Pd-Tsph code was written during the 
' 90s by C. Chiosi, G. Carraro. C. L ia and C. Dalla Vecchia 
tearraro et Il998l : iLi a et al.1 l2002h . Over the years, many 
researchers added their contribution to the development of 
the code. In its original release, the Pd-Tsph was a basic 
Tree + SPH code, written in Fortran 9 Q, co nceptually simi- 
lar to Tree-SPH bv lHernauist Katd (| 19891 ). Schematica lly, 
Pd-Tsph used an early formulation of SPH teen j \l99d ) to 
solve the equa t ions o f motion for the gas component, and the 
iBarnes HutI (| 19861 ) Tree algorithm to compute the gravita- 
tional interactions. 



boundary conditions, and the parallelization of the code. 
Sect. [5] contains and discusses numerous tests aimed at 
assessing the performance of the code. Sect l6] summarizes 
some concluding remarks. Finally, Appendices A and B 
contains some technical details related the N-dimensional 
spline kernel (Appendix A) and the equations of motion 
and energy conservation (Appendix B). 

2. Description of the code: Gravity 

Gravity is the leading force behind the formation of cos- 
mic structures, on many scales of interest, from galaxy 
clusters down to individual stars and planets. Classical 
gravitation is a well understood interaction. As long as 
General Relativity, Particle Physics or exotic treatments 
such as Modified Dynamics are not considered 0, it is de- 
scribed by Newton's law of universal gravitation: 

F.,=G^^^{r,-n), (1) 

where G is the gravitational constant. Given a density field 
p(r), the gravitational potential $ is obtained via Poisson 
equation: 

V^^{r) =47rGp(r), (2) 

and Fij = V<l>(r). 

The gravitational force exerted on a given body (i.e., 
particle) by the whole system of bodies within a simula- 
tion can be obtained by the vectorial summation of all the 
particles' contributions (without considering, for the mo- 
ment, the action of the infinite region external to the com- 
putational volume; this can indeed be taken into account 
using periodic boundary conditions, see Sect. 14. 3j) . This is 
simply the straightforward application of Eq. [TJ Anyway, 
in practice this approach is not efficient, and may also 
lead to artificial divergences, as explained in the following 
Sections. 

2.1. Softening of the gravitational force 

Close encounters between particles can cause numerical 
errors, essentially because of the time and mass resolu- 
tion limits. Moreover, when dealing with continuous fiu- 
ids rather than with single, isolated objects, one tries to 
solve Eq. [2] rather than Eq. [H and must therefore try to 
model a smooth density field given a distribution of par- 
ticles with mass. In addition, dense clumps of particles 
may also steal large amounts of computational time. To 
cope with all these issues, it is common practice to soften 
the gravitational force between close pairs of bodies: if 
the distance between two particles becomes smaller than 
a suitable softening length e, the force exerted on each 
body is corrected and progressively reduced to zero with 
decreasing distance. 

^ A relativistic formulation of both gravitational and hydro- 
dyn amical interact ions is possible (for a general summary see 
e.g. lRosswodl2009h , but is generally unessential in problems of 
cosmological structure formation. 
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Different forms of the force softening function can be 
used. A possible expression is given by the following for- 
mula: 



F{r,e) 
Grrii 
'l/e2r4 

1/lr 



(3) 



3« - f 



1 

15n2 



if < < 1, 
if 1 < < 2, 
if > 2, 



where r is the distance between particles i and j, and u = 
\r\/e. This expression for the force softening corresponds 
(via the Poisson equation) t o a densit y distribu tion k ernel 
function proposed by M onaghan fc Lattanziol ( 1985[ ) and 



widely adopted in Smoothed Particles Hydrodynamics al- 
gorithms (see Sect. [3]) H: 



W{r, e) = —- X < 



i(2- 



if < u < 1, 
if 1 < u < 2, 
if u > 2. 



(4) 



Note that in Eq. 2] the softening length e is assumed to 
be the same for the two particles i and j. In the more 
general situation in which each particle carries its own e, a 
symmetrisation is needed to ensure energy and momentum 
conservation: this can be achieved either using e = elj = 
(ei + ej)/2, or by symmetrizing the softened force after 
computing it with the two different values ei and Cj. 

The softening lengths can be fixed in space and time, 
or may be let vary with local conditions. If the soften- 
ing length is kept constant, the choice of its numerical 
value is of primary importance: for too small a soften- 
ing length it will result in noisy force estimates, while 
for too large a value it will systematically bias the force 
in an unp hysical manne r (Merritt 1996|;_ 



7T"--^^ ^. 1 9961: IRQme o| JT998: 

Athanassoula et al I I2QQQI : [Price MonaghanI 120071 . and 



see also Sect. 15. 9p . Unfortunately, the "optimal" value for 
the softening length depends on the particular system be- 
ing simulated, and in general it is not known a priori. A 
standard solution is to assign to each particle a softening 
length proportional to its mass, keeping it fixed through- 
out the whole simulatiorjfl, or letting it vary with time, 
or redshift, if a cosmological background is included. A 
clear advantage of keeping e fixed in space is that energy 
is naturally conserved, but on the other hand the smallest 



^ Throughout the paper, we only use the 3-dimensional for- 
mulation of kernels. See Appendix A for the 1-D and 2-D forms. 

^ Recently, Shirokov (2007) pointed out that since particles 
are of unequal mass, and hence unequal softening lengths, one 
should actually compute the pairwise gravitational force and 
potential by solving a double integral over the particle volumes. 
Therefore, the computation of the interactions using the classic 
approach is likely not accurate. Given the practical difficulty 
in evaluating those integrals, they also provide a numerical 
approximation. 



resolvable length scale is fixed at the beginning of the sim- 
ulation and remains the same in the whole spatial domain 
independently of the real evolution of the density field. 
A collapsing or expanding body may quickly reach a size 
where the flow is dominated by the softening in one re- 
gion, while in another the flow may become unphysically 
point-like. 

Obviously, the accuracy can be greatly improved if e 
is let var y according; to the local particle number density 
(see e.g. IPehnenl 200l|). Moreover, in principle, if parti- 
cles in a simulation are used to sample a continuous fluid 
(whose physics is determined by the Navier- Stokes or the 
Boltzmann equations) the properties of such points should 
always change accordingly to the local properties of the 
fluid they are sampling, to optimize the self-adaptive na- 
ture of Lagrangian methods. On the other hand, if parti- 
cles represent discrete objects (single stars, or galaxies in 
a cluster, etc.), their softening lengths might perhaps be 
considered an intrinsic property of such objects and may 
be kept constant, depending on their mass and not on 
the local density of particles. In cosmological and galaxy- 
sized simulations, gas and Dark Matter particles are point- 
masses sampling the underlying density field, and stellar 
particles represent clusters of thousands of stars prone to 
the gravitational action of the nearby distribution of mat- 
ter; thus a fully adaptive approach seems to be adequate 
to describe the evolution of these fluids. However, it can 
be easily shown that simply letting e change freely would 
result in a poor conservation of global energy and momen- 
tu m, even if in some cases the errors could be sma ll (see 
e.g. IPrice fc M onaghanl l2007l : IWetzstei n eralll2008h . 

To cope with this, EvoL allows for the adaptive evo- 
lution of individual softening lengths, but includes in the 
equations of motion of particles self-consistent correct- 
ing additional terms. S uch terms can be derived (see 



Price fc Monaghanll2007l ) starting from an analysis of the 
Lagrangian describing a self-gravitating, softened, colli- 
sionless system of N particles, which is given b}{l 



N 



(5) 



where ^ is the gravitational potential of the i-th particle, 

N 



^{Vi) = -G^mj(l){\rij\,ei), 
and (/) is a softening kernel (r^ 



(6) 



Vi — Vj). Assuming 
a density distribution described by the standard spline 
kernel Eq. [3J the latter becomes 



^ A clear advantage of using the Lagrangian to derive the 
equations of motion is that, provided the Lagrangian is appro- 
priately symmetrised, momentum and energy conservation are 
guaranteed. Note that this derivation closely matches the one 
described in Sect.[3]to derive the variational formulation of the 
SPH formalism and the so-called V/i correcting terms. 
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0(r,e) 



-\I\t\ 



10^ 



10^ 



10^ 



30^ 



(7) 

if < < 1, 

I + -r^l if 1< < 2, 

5 15n-i — ' 

if > 2, 



Approximating the integral to a summation over par- 
ticles for practice purposeqH, one obtains 



where again u = \r\/e. Note that the force softening, Eq. 
m is obtained taking the spatial derivative of Eq. [HI 

The equations of motion are obtained from the Euler- 
Lagrange equations, 



d 
dt 



dL 
dvi 



dL_ 



0, 



which give 

dVi dL 
dt dVi ' 

The gravitational part of the Lagrangian is 



m 



N 

j k 



(8) 



(9) 



(10) 



where for the sake of clarity (e^) = ^(Ir^j |, e^); swapping 
indices, the partial derivative dLgrav/dvi results 



(14) 



where i marks the "active" particle that lays at the posi- 
tion at which the function is being evaluated, and j are its 
neighbouring particles (or simply neighbours). The error 
by doing so de pends on the diso rder of particles and is in 
general 0{h'^) (|Monaghanlll992h . Thus, 



ni = Y^W{\rijlei), 

j 



(15) 



where it has been put A = n; so, there is no need to 
know in advance the value of the density of neighbouring 
particles to compute the density of the active one. 
To link e and n, a safe prescription is to use 



V I — 

. rii 



1/3 



(16) 



where 7^ is a dimensionless parameter. With this law, the 
weighted number of particles within a softening sphere is 
tentatively held constant, i.e. 



dL 



grav 



dvi 
G 



(11) -{2e,fni=Nn 



d^jkicj) d\rjk\ d^jk{ej) dej 



d\rjk\ 



where 

d\rjk\ ^ Vj - Vk 
dvi h 



rk 



{Sji — Ski). 



dvi 



(12) 



To obtain the required terms in the second addend on 
the right part of Eq. [121 the softening length must be re- 
lated to the particle coordinates so that e = e(n), where n 
is the number density of particles at particle i location. To 
this aim, one can start from the general interpolation rule 
that any function of coordinates A{r) can be expressed in 
terms of its values at a disordered set of points, and the 
integral interpolating the function A{r) can be defined by 

Ai^t{r) = J A{r')W{r - r', e)dr', (13) 

where W is the density kernel Eq. [H This is the same rule 
at the basis of the SPH scheme (see Sect. [3]) El- 

^ Note that in principle W should be defined over the whole 
space, as in the first SPH calculations bv lGingold MonaghanI 
(|l977l l where a Gaussian-shaped kernel was adopted. For prac- 
tical purposes, anyway, its domain is made compact putting 
VK = for distances greater than r/e, with 77 > 0, from the 
position of the active particle. 



(17) 



with Nnei,id^7r{2r])^ . Numerical experiments have shown 
that a safe ch oice is 1.2 < 77 < 1.5, corresponding to 60 < 
Nnei,id < 110 (|Price fc MonaghanlBoO?! 
The term needed in Eq. [12] is 



dcj dcj duj 



dvi 



drij dvi 



(18) 



The first factor is given by the derivative of Eq. [161 



dej_ 
drin 



3nn 



(19) 



whereas the second factor is the spatial derivative of Eq. 
[T5l which results: 



drij 1 



dWj,{e,) 
dvi 



(20) 



^ IPrice fc MonaghanI (|2007l ^ use mass the density in place 
of number density to achieve the desired solutions. The for- 
mulation in terms of number density is equivalent and can be 
advantageous when dealing with particles of different masses. 
Also note that this is a reformulation of the well known SPH 
density summation, see Sect. [3] 
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where it was defined 



1 



dei ^ dWij{ei) 



dn. 



Rearranging and putting G . 

dLnrav 



1, one finally finds 



(21) 



(22) 



-m. 



2 

dWij{ei 
T,: dvi 



ij dWijici 



T 



drj 



where 



d(j)ij{ei) 



(23) 



The d(j)/de terms used in Eq.[23]can be tabulated together 
with the other kernel functions: 



proposed an iterative method in which, beginning from a 
starting guess for both n and e, the solution of the equa- 
tion 



/(e) = \ni{e,) 



(25) 



where nsum,i{^i) is the mass density computed via sum- 
mation over neighbouring particles (Eq. [T5|) and ni{ei) is 
the density obtained from Eq. [161 is searched by means 
of an iterative procedure, that can be solved adopting a 
classical bisection scheme, or more efficient routines such 
as the Newton-Raphson method. In this case, a loop is 
iterated until \ei^riew - ^i\/^i,init < Itol, where ^tol is 
a tolerance parameter ~ 10~^ — 10~^, ei^init is the value 
of e for the particle i at the beginning of the procedure, 
Ci is its current value, and Ci^new = — f{^i)/f'{^i)] the 
derivative of Eq. [25] is given by 



j 



(26) 



de 



(24) 



l/e2[ 




-2^2. 



if < < 1, 
if 1 < < 2, 
if ^1 > 2, 



where, as usual, u = \r\/e. Finally, the dW/dr and dW/de 
terms are easily obtained deriving the explicit expression 
ofW{\r\,e). 

In Eq. 22 (and Eq. 12), the first term is the classi- 
cal gravitational interaction, whereas the second term is 
a new, extra-term which restores energy conservation for 
varying softening lengths. Since ^ is defined as a neg- 
ative quantity for positive kernels, this term acts as a 
negative pressure, increasing the gravitational force. To 
obtain the ^ and T correcting terms, each particle i 
must perform a loop over the other particles, summing 
the contributions from the ones that satisfy the criterion 
|ri_^-| < 2 X max(ei,e_;-)EI. 

The relation between e and n, defined by Eq. [TH 
leads to a non-linear equation which can be solv e d self - 
consistently for each particle. Price fc MonaghanI ( 2007 ) 



^ Thus, all particles, and not only SPH particles, need to find 
their neighbours to compute these terms (see Sect. 12.2]) . While 
the Q term described in Sect. [3] is only needed for the SPH 
algorithm, and therefore only SPH particles concur to com- 
pute it (considering only SPH neighbours in turn), the ^ and 
T gravitational terms are required for any kind of gravitating 
particle, and must be computed looping on all neighbouring 
particles, both SPH and not. The time lost in the neighbour 
searching routine for non-gaseous particles can be at least par- 
tially balanced adopting the individual time-stepping scheme 
(see Sect. KTB 



In this case, large softening lengths in low- 
density regions result in larger individual time-steps for parti- 
cles belonging to those regions. 



To increase the efficiency of this process, a predicted 
value for both e and n can be obtained at the beginning 
of each time-step using a discretized formulation of the 
Lagrangian continuity equation. 



dn ,^ 
- = -n(V..). 



(27) 



Such formulation can be obtained taking the time deriva- 
tive of Eq. [151 which results 



dni 

~dt 

J 

while 

dci dci dni 
dt dni dt 



^^{vi-Vj)-ViW{rij,e), 



(28) 



(29) 



Combining Eq. [28] with Eq. [271 one can also see that 
the velocity divergence at i particle location is given by 



V-Vi 



vj)-ViWirij,e). 



(30) 



The adoption of this adaptive softening length formal- 
ism, with the correcting extra-terms in the equation of 
motion, results in small errors, always comparable in 
magnitude to the ones f ound with the "optimal" e (see 
Price fc Monaghan|[2QQ7l .their Figs. 2, 3 and 4). 

At the beginning of a simulation, the user can select 
whether he/she prefers to adopt the constant or the adap- 
tive softening lengths formalism, switching on or off a suit- 
able flag. 
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2.2. Hierarchical oct-tree structure 

A direct summation of the contribution of all particles 
should in principle be performed for each particle at each 
time-step to correctly compute the gravitational interac- 
tion, leading to a computational cost increasing with N'^ 
{N being the number of particles). A convenient alter- 
native to this unpractical approach are the so-called tree 
structures, in which particles are arranged in a hierar- 
chy of groups or "cells" . In this way, when the force on a 
particle is computed, the contribution by distant groups 
of bodies can be approximated by their lowest multipole 
moments, reducing the computational cost to evaluate the 
total force to 0{NlogN). In the classical iBarnes fc Hut 
( 19861 ) scheme, the computational spatial domain is hi- 
erarchically partitioned into a sequence of cubes, where 
each cube contains eight siblings, each with half the side- 
length of the parent cube. These cubes form the nodes of 
an oct-tree structure. The tree is constructed such that 
each node (cube) contains either a desired number of par- 
ticles (usually one, or one per particles type - Dark Matter, 
gas, stars), or is progenitor of further nodes, in which case 
the node carries the monopole and quadrupole moments 
of all the particles that lie inside the cube. The force com- 
putation then proceeds by walking along the tree, and 
summing up the appropriate force contributions from tree 
nodes. In the standard tree walk, the global gravitational 
force exerted by a node of linear size / on a particle i is 
considered only \i r > l/O^ where r is the distance of the 
node from the active particle and is an accuracy param- 
eter (usually < 1): if a node fulfills the criterion, the tree 
walk along this branch can be terminated, otherwise it is 
"opened" , and the walk is continued with all its siblings. 
The contribution from individual particles is thus consid- 
ered only when they are sufficiently close to the particle 
i. 

To cope with som e problems pointed out by 
Salmon fc Warrenl ( 1994 ) about the worst-case behaviour 
of this standard criterion for commonly employed opening 
angles, Dubinski (1996.) introduced the simple modifica- 
tion 



r>l/e^ 5, 



(31) 



where the 6 is the distance of the geometric center of a 
cell to its center of mass. This provides protection against 
pathological cases where the center of mass lies close to 
an edge of a cell. 

In practice, this whole scheme only includes the 
monopole moment of the force exerted by distant groups 
of particles. Higher orders can be included, and the 
common practice is to include the quadrupole mo- 
ment corrections (EvoL indeed includes them). Still 
higher multipole orders would result in a worst perfor- 
manc e without significant ^ains in computational accu- 
racy (jMcMillan fc Aarseth|[l993[ ). 

Note that if one has the total mass, the center of mass, 
and the weighted average softening length of a node of the 
oct-tree structure, the softened expression of the gravita- 



tional force can be straightforwardly computed treating 
the cell as a single particle if the opening criterion is ful- 
filled but the cell is still sufficiently near for the force to 

need to be softened. 

As ffrst suggested by iHernquist fc Katz ( 19891 ). the 
tree structure can be also used to obtain the individual 
lists of neighbours for each body. At each time step each 
(active) particle can build its list of interacting neighbour- 
ing particles while walking the tree and opening only suf- 
ffciently nearby cells, until nearby bodies are reached and 
linked. 



3. Description of the code: Hydrodynamics 

Astrophysical gaseous plasmas can generally be reason- 
ably approximated to highly compressible, unviscous ffu- 
ids, in which anyway a fundamental role is played by 
violent phenomena such as strong shocks, high energy 
point-like explosions and/or supersonic turbulent motions. 
EvoL follows the basic gas physics by m eans of the 



Smoothed Parti cles Hydrodynamics (SPH, ILucvI 



Monaghanlll992r). i n a m odern formulation based on the 
review bv iRosswogl ( 2QQ9[ ). to which the reader is referred 
for details. Anyway, some different features are present in 
our implementation and we summarize them in the follow- 
ing, along with a short review of the whole SPH algorithm. 

To model ffuid hydrodynamics by means of a discrete 
description of a continuous the ffuid, in the SPH approach 
the properties of single particles are smoothed in real space 
through the kernel function VK, and thus weighted by the 
contributions of neighbouring particles. In this way, the 
physical properties of each point in real space can be ob- 
tained by the summation over particles of their individual, 
discrete properties. Note that, on the contrary of what 
happens when softening the gravitational force (which is 
gradually reduced to zero within the softening sphere), in 
this case the smoothing sphere is the only "active" region, 
and only particles inside this region contribute to the com- 
putation of the local properties of the central particle. 

Starting again from the interpolation rule described in 
Sect. 12. H but replacing the softening length e by a suitable 
smoothing length the integral interpolating the function 
A{r) becomes 



Ai{r) = J A{r')W{r - r^ h)dr', 



(32) 



(the kernel function W can be the density kernel Eq. [4]). 
The relative discrete approximation becomes 



Pj 



(33) 



and the physical density of any SPH particle can be com- 
puted as 



Pi 



Y,m,W{\r\M) 



(34) 
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(here and throughout this Section, it is imphed that ah 
summations and computations are extended on SPH par- 
ticles only). 

A differentiable interpolating function can be con- 
structed from its values at the interpolation points, i.e. 
the particles: 



VA(r) = ^mjVW{r - r^ h). 



(35) 



Anyway, better accuracy is found rewriting the formu- 
lae with the density inside operators and using the general 
rule 

p\/A = V(M) - A\/p (36) 

(e.g., 'Monaghanl ll992[ ). Thus the divergence of velocity is 
customarily estimated by means of 

\/-Vi = ^mjlivj - v^) ■ ViWij]/pi, (37) 

3 

where ViWij denotes the gradient of W{r — r' ^h) with 
respect to the coordinates of a particle i. 

The dynamical evolution of a continuous fluid is gov- 
erned by the well known laws of conservation: the con- 
tinuity equation which ensures conservation of mass, the 
Euler equation which represents the conservation of mo- 
mentum, and the equation of energy conservation (plus 
a suitable equation of state to close the system). These 
equations must be written in discretized form to be used 
with the SPH formalism, and can be obtained by means of 
the Lagrangian analysis, following the same strategy used 
to obtain the gravitational accelerations. 

The continuity equation is generally replaced by Eq. 
£|. Like in the case of the gravitational softening length 
e, the smoothing length h may in principle be a con- 
stant parameter, but the efficiency and accuracy of the 
SPH method is greatly improved adapting the resolution 
lengths to the local density of particles. A self-consistent 
method is to adopt the same algorithm described in Sect. 
2.1\ i.e. obtaining h from 



dh dh dn 
dt dn dt ' 



(39) 



where of course n is the local number density of SPH par- 
tides only^ and relating to n by requiring that a fixed 
number of kernel- weighted particles is contained within a 
smoothing sphere, i.e. 



hi 



1/3 



Alternatively, it could be written as 



dt 



(40) 



(38) 



which is advantageous in case the simulated fluid has bound- 
aries or edges, but has the disadvantage of not conserving the 
mass exactly. 



Note that in this case, while still obtaining the mass den- 
sity by summing Eq. [34j one can compute the number 
density of particles at particle i's location, i.e.: 



E 

3 



(41) 



which is clearly formally equivalent to Eq. [rfpl . 

It is worth recalling here that in the first versions of 
SPH with spatially-varying resolution, the spatial varia- 
tion of the smoothing length was generally not considered 
in the the equations of motion, and this resulted in secu- 
lar errors in the conservation of entropy. The inclusion of 
extra-correcting terms (the so-called Vh terms) is there- 
fore important to ensure the conservat i on of both energy 
and entropy. For example. ISerna] et al.l (|l996l ) studied the 
pressure-driven expansion of a gaseous sphere, finding that 
while the energy conservation is generally good even with- 
out the inclusion of the corrections (and sometimes get 
slightly worse if these are included!), errors up to ~ 5% can 
be found in the entropy conservation, while all this does 
not occur with the inclusion of the extra-terms. Although 
the effects of entropy violation in SPH codes are not com- 
pletely clear and need to be analysed in much more de- 
tail, especially in simulations where galaxies are formed in 
a cosmological framework, there are many ev idences that 
the pr oblem must be taken into consideration. [Alimi et al] 
(|2QQ3l ) have analysed this issue in the case of the collapse of 
isolated objects and have found that, if correcting terms 
are neglected, the density peaks associated with central 
cores or shock fronts are overestimated at a :^ 30% level. 

These Vh terms were first introduced explicitly 
([Nelson fc Papaloizoul Il994h . with a double-summation 
added to the canonical SPH equ ations. Lat er, an 



implicit formulation was obtained (|Monaghanl l2QQ2l : 



Springel fc Hernquistll2QQ2l ). starting from the Lagrangian 



derivation of the equations of motion and self-consistently 
obtaining correcting terms which accounts for the varia- 
tion of h. Obviously such terms are formally similar to 
those obtained for the locally varying gravitational soft- 
ening lengths (Sect. 12.11) . 

Following iMonag hinl (|2002l ). the Lagrangian for non- 
relativistic fluid dynamics can be written as 

'1 



L 



P 



u dV, 



(42) 



(the gravitational part is not included here since SPH does 
not account for gravity), which in the SPH formalism be- 
comes 



-'SPH 



1 2 



(43) 



Some authors (e.g. IHu k Adam j l20Q5l : lOtt k Schn^tte3 
l2003l ) proposed a different "number density" formulation of 
SPH in which pi — miXm, and Ui is obtained via Eg. 1411 While 
this can help in resolving sharp discontinuities and taking into 
account the presence of multi-phase flows, it may as well lead to 
potentially disastrous results if mixed unequal mass particle are 
not intended to model density discontinuities but are instead 
used as mass tracers in a homogeneous field. 
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As already made for Eq. [8l the equations of motion of 
a particle can be obtained from the Euler-Lagrange equa- 
tions, giving 



~dt 



E 



du\ dpj 



dp J dvi 



(44) 



To obtain the term dpj/dvi^ one can note that using the 
density p given by Eq.[34land letting the smoothing length 
vary according to Eq. [40] one gets 



d 1 

^ = ^ ^ mkViWik{hi)5ij - rriiVjWijihj), 



(45) 



where 6ij is the Kronecker delta function, \/iW is the gra- 
dient of the kernel function W taken with respect to the 
coordinates of particle i keeping h constant, and 



k ^ 



(46) 



is an extra-term that accounts for the dependence of the 
kernel function on the smoothing length. Note that Q is 
formally identical to the T term introduced in Sect. [21 
replacing e with h. This is obvious since the underlying 
mathematics is exactly the same; both terms arise when 
the motion equations are derived from the Lagrangian for- 
mulation. 

In case of particles with different mass, the term 
dpj/dvi is 



k 



Cj/rrik 



where 
dhi 



and 



drij 



E 



dhi 



-Sik), (47) 



(48) 



n* = 1 



dhi 



^aWy(/ii) 

j 



dhi 



(49) 



(Price 2009, private communication). These terms can be 
easily computed for each particle while looping over neigh- 
bours to obtain the density from Eq. [Ml 

The term du/dp in Eq. [Hj can be derived from the 
first law of thermodynamics, dU = dQ — PdV ^ dropping 
the dQ term since only adiabatic processes are considered. 
Rewriting in terms of specific quantities, the volume V 
becomes volume "per mass", i.e. and dV = d{l/p) = 
-dp/p^. Thus, 



du 



—dp, 



(50) 



and, if no dissipation is present, the specific entropy is 
constant, so 



dt ) 



P 
7^' 



(51) 



Inserting Eqs. [511 and [TTl into Eq. [33J the equations of 
motion finally become 



dvi ^ 
— — = > rrij X 
dt ^ ^ 



(52) 




The equation for the thermal energy equation is 



dui 
~di 



^2 



Ci/rrij 



(53) 

Vi)-\/iWij{hi) . 



It can be shown that if all SPH particles have the same 
mass, Eqs. [53l and [54l reduce to the standard Eqs. 2.10 and 
2.23 in Monaghan (2002)0- 

Equations HH (density summation) and Eq. [40l form a 
non-linear system to be solved for h and n, adopting the 
same method described in Sect. 12.11 

A drawback of this scheme can be encountered when 
a sink of internal energy such as radiative cooling is in- 
cluded. In this case, very dense clumps of cold particles 
can form and remain stable for long periods of time. Since 
the contribution of neighbouring particles to the density 
summation (Eq. [41]) is weighted by their distance from 
the active particle, if such clump is in the outskirts of the 
smoothing region then each body will give a very small 
contribution to the summation, and this will result in an 
unacceptably long neighbour list. In this cases, the adop- 
tion of a less accurate but faster method is appropriate. 
The easiest way to select h is to require that a constant 
number of neighbouring particles is contained within a 
sphere of radius rjh (perhaps allowing for small devia- 
tions). This type of procedure was generally adopted in 
the first SPH codes, and provided sufficiently accurate re- 
sults. If the scheme in question is adopted (EvoL may 



As noted by [Schuessler Schmitt[ (|l98l[ l. the gradient of 
the spline kernel Eq. [4] can lead to unph ysical clustering of 
particles. To prevent this, Monaghan (2000) introduced a small 
artificial repulsive pressure-like force between particles in close 
pairs. In EvoL a modifie d form of the kernel grad i ent is instead 
adopted, as suggested bv [Thomas Couchman] (|l992[ ): 



dW 
dr 



(4 if0<i^<2/3, 
3u(4-3u) if2/3<u< 1, 



47r 



3(2- 
10 



if 1 < 1^ < 2, 
if 1^ > 2, 



(54) 



with the usual meaning of the symbols. 



E. Merlin et al.: EvoL I 



9 



choose between the two algorithms), the terms \/h are ne- 
glected, because the relation Eq. [40] is no longer strictly 
valid. In the following, we will refer to this scheme as to 
the standard SPH scheme^ and to previous one based on 
the n — h (or p — h relation) as to the Lagrangian SPH 
scheme. 

A test to check the ability of the different algorithms 
to capture the correct description of an (ideally) homo- 
geneous density field has been performed, using parti- 
cles of different masses mixed together in the domain 
< X < 1,0 < y < 0.1,0 < z < 0.1, to obtain the 
physical densities p = 1 if x < 0.5 and p = 0.25 if x > 0.5. 
To this aim, a first set of particles were displaced on a 
regular grid and were given the following masses: 

m = 0.9 X 10"^ for x < 0.5; m = 2.25 x 10"^ for x > 0.5. 

Then, a second set of particles were displaced on another 
superposed regular grid, shifted along all three dimensions 
with respect to the previous one by a factor Ax = Ay = 
Az = 5 X 10~^. These particles were then assigned the 
following masses 

m = 0.1 X 10"^ for x < 0.5; m = 0.25 x 10"^ for x > 0.5. 

Then, four different schemes were used to compute the 
density of particles: 

— standard SPH, mass density summation 

— standard SPH, number density summation 

— Lagrangian SPH, mass density summation 

— Lagrangian SPH, number density summation 

where the expressions "mass density" and "number den- 
sity" refer to the different scheme adopted to compute the 
density of a particle, i.e. to Eq. [SH and [HJ respectively. 
Looking at Fig. [H it is clear that the discrete nature of 
the regular displacement of particles gives different results 
depending on the adopted algorithm. The mass density 
formulation is not able to properly describe the uniform 
density field: low-mass particles strongly underestimate 
their density, on both sides of the density discontinuity. 
The situation is much improved when a number density 
formulation is adopted. Little differences can be noted be- 
tween the determinations using the constant neighbours 
scheme and the density-/^ relation. 

3.1. On symmetrisation 

To ensure the conservation of energy and momenta, a 
proper symmetrisation of the equations of motion is re- 
quired. It is worth noticing that while in the first formula- 
tions of SPH such symmetrisation had to be imposed "ad 
hoc" , in the Lagrangian derivation it is naturally built-in 
without any further assumption. 

Anyway, the symmetrisation is strictly necessary only 
when pairwise forces act on couples of particles. Thus, it is 
not needed when obtaining the density from Eq.[4T]and[34j 
or when the internal energy of a particle is evolved accord- 
ing to Eq. [SH Indeed, symmetrisation of the energy equa- 
tion may lead to unphysical negative temperatures. For 



0.5 




0.5 




0.5 



0.5 



0.5 

X 



0.5 

X 



Fig. 1. Determination of the Initial density adopting four 
different SPH algorithms. Left to right, top to bottom: 
standard SPH, mass density summation; Lagrangian SPH, 
mass density summation; standard SPH, number density 
summation; Lagrangian SPH, number density determina- 
tion. Blue dots: high mass particles; red crosses: low mass 
particles. 



example when cold particles are pushed by a strong source 
of pressure: in this case, the mechanical expansion leads 
to a loss of internal energy which, if equally subdivided 
between the cold and hot phase, ma y exceed the small en- 



Springel fc Hernquist 



ergy budget of cold bodies (see e.g 
20021) . 

Since symmetrisation is not needed in the density sum- 
mation algorithm, the smoothing length of neighbouring 
particles is not required to compute the density of a par- 
ticle i. This allows an exact calculation simply gathering 
position and masses of neighbours, which are known at 
the beginning of the time-step. On the other hand, af- 
ter the density loop each particle has its new smooth- 
ing length h and it is now possible to use these values 
in the acceleration calculations, where symmetrisation is 
needed. Anyway, while in the previous case only parti- 
cles within 2hi were considered neighbours of the particle 
z, in this latter case the symmetrisation scheme requires 
a particle j to be considered neighbour of particle i if 



< 2 X max(/i^, hj) 



12 



Note that with this scheme two routines of neighbour 
searching are necessary: the first one in the density sum- 
mation, and the second in the force calculation. In prac- 
tice, when the evaluation of density i s being performed, 
a gather view of SPH is adopted (see iHernquist fc Katz 



It would also be possible to compute the interactions by 
only considering neighbours within hi (or a). To do so, each 
particle should "give back" to its neighbours its contribution to 
the acceleration. Anyway, this approach is a more complicated 
with individual timesteps (where one has to compute inactive 
contributions to active particles), and in a parallel architecture. 
We therefore decided to adopt the fully symmetric scheme in 
this release of EvoL 
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Il989h : thus, each particle i searches its neighbours using 
only its own smoothing length /i^, until the convergence 
process is completed. In this way, if a particle finds too 
high or too low a number of neighbours (and therefore 
must change its smoothing length to repeat the search) 
it is not necessary to inform the other particles (includ- 
ing those belonging to other processors) about the change 
that is being made. This would in principle be necessary 
if a standard gather /scatter mode is adopted. Thus, each 
particle evaluates its density independently, adjusting its 
smoothing length on its own, and no further communica- 
tion among processors is necessary, after the initial gath- 
ering of particles positions. Moreover, the search is ini- 
tially performed using an effective searching radius which 
is slightly larger than the correct value 2hi. In this way, 
the neighbour list will include slightly more particles than 
effectively needed. However, if at the end of the loop, the 
smoothing length has to be increased by a small factor, it 
is not be necessary to repeat the tree- walk to reconstruct 
the neighbour list. Of course, if the adaptive softening 
length scheme is adopted the tree- walk will be performed 
to upgrade h and p for SPH particles only, and e and ptot 
for all particles. 

Then, when forces (hydrodynamical forces and, if the 
adaptive softening length scheme is adopted, gravitational 
forces too) are evaluated, a new search for neighbours must 
be made, in which a particle is considered neighbour to 
another one if the distance between the two is lower than 
the maximum between the two searching radii. During this 
second tree-walk, the gravitational force is evaluated to- 
gether with the construction of the particle's neighbour 
list. Finally, correcting terms are computed for hydro- 
dynamics equations of motion {Vh terms) and, if neces- 
sary, for gravitational interactions (if softening lengths are 
adaptive), summing the neighbouring particles' contribu- 
tions. 



3.2. Smoothing and softening lengths for SPH particles 



SPH particles under self-gravity have two parameters 
characterizing their physical properties: the softening 
length e, which tunes their gravitational interactions with 
nearby particles, and the smoothing length /i, used to 
smooth their hydrodynamical features. 

In principle there is no need for these two quanti- 
ties to be related to one another, because they refer to 
two distinct actions, and have s omehow opposit e func- 
tions as noted above. However, iBate BurkertI (|l997h 
claimed that a quasi-stable clump of gas could became 
un-physically stable against collapse if e > /i, because in 
this case pressure forces would dominate over the strongly 
softened gravity, or on the other hand, if e < /i, col- 
lapse on scales smaller than the formal resolution of SPH 
may be induced, causing artificial fragmentation. They 
recommend that gravitational and hydrodynamical forces 
should be resolved at the same scale length, imposing e = 
h (requiring both of them to fixed, or e to be adaptive like 



/i, with the introduction of the correcti ng terms described 
above). In other studies (e.g. Thack er fc Couchman 20001) 
the softening length is fixed, but the smoothing length is 
allowed to vary imposing a lower limit close to the value of 
the softening length (this is a typical assumption in SPH 
simulations of galaxy format i on). 

Anyway, I Williams et al. ([2004') have studied the ef- 
fects of such procedure, finding that in many hydrody- 
namical studies this can cause a number of problems. The 
adoption of large values of h can result in over-smoothed 
density fields, i.e. a decrease in the computed values of 
density p and of density gradients Vp, which in turn un- 
physically decreases the computed pressure accelerations 
in the Euler equation. Problems with angular momentum 
transfer may also arise. Therefore, they suggest to allow 
h to freely vary. This eventually yields improved accuracy 
and also, contrary to expectations, a significant saving in 
computational time (because the number of neighbours is 
kept fixed instead of largely increasing as it may happen in 
dense environments with fixed h). Moreo ver, they pointed 
out that, in contrast with the claims by Bate fc BurkertI 



(|199Z ). the smoothing length should not be kept equal 
to the softening length, because in many physical situa- 
tions (shocks for example) hydrodynamical forces domi- 
nate over gravity and likely need to be properly resolved 
on size scales smaller than the softening length. 

Finally, a numerical problem arises when trying to 
keep e = h in cosmological simulations. The presence of 
a Dark Matter component makes it impossible to keep 
constant the number of both hydrodynamical and gravi- 
tational neighbours. For example, a collapsed Dark Matter 
halo could retain few gaseous particles (because of shocks 
and/or stellar feedback), and these will thus have many 
gravitational neighbours but few hydrodynamic neigh- 
bours with a unique softening/smoothing length. The op- 
posite problem could arise in the case of very dense clumps 
of cold gas that may form behind shocks, without a cor- 
responding concentration of Dark Matter. 

Thus, in EvoL both e and h are let free to vary in- 
dependently from one another. Future tests may be done 
trying to keep the two parameters linked, for example al- 
lowing for some variations in both their ratio e/h (impos- 
ing it to be around^ but not exactly, 1) and in the number 
of hydrodynamical and gravitational neighbours. 

3.3. Discontinuities 
3.3.1. Artificial viscosity 

The SPH formalism in its straight form is not able to 
handle discontinuities. By construction, SPH smooths out 
local properties on a spatial scale of the order of the 
smoothing length h. Anyway, strong shocks are of primary 
importance in a number of astrophysical systems. Thus, 
to model a real discontinuity in the density and/or pres- 
sure field ad hoc dissipational, pressure-like terms must 
be added to the perfect fiuid equations, to spread the dis- 
continuities over the numerically resolvable length. This 
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is usually achieved by introducing an artificial viscosity^ 
a numerical artifact that is not meant to mimic physi- 
cal viscosity, but to reproduce on large scale the action 
of smaller, unresolved scale physics, with a formulation 
consistent with the Navier-Stockes terms. 

Among the many formulations that have been pro- 
posed, the most commonly used is obtained adding a term 
Hij to the Eqs. [53land[5^ so that 



dt ~ ^ 



rrij X 



(55) 




ViWij{hi)- 



ViWij{hj)+IlijViWi, 



and 

dui 
'dt 



(56) 



E 



rrij 



El 

^2 



1 + ^) {v,-v,)-V,W.,{h,) 



-IiijWij\ViWij\ 



Here Wij (vij • rij)/\rijl 



Hi 



if Vij ■ Tij < 0, 



if Vi 



>0, 



hvi 



WiWij{hi)(l 



+ViWijihj) 




(57) 

(58) 

(59) 
(60) 



(lMonaghanlll988l ). Cij is the average sound speed, which 



is computed for each particle from a suitable equation 
of state (EoS), like the ideal monatomic gas EoS P = 
(7 — l)pu, using Ci = Y^7(7 — l)ui (7 is the adiabatic 
index); the parameter r] is introduced to avoid numeri- 
cal divergences, and it is usually taken to be r] = O.lh. 
a and /3 are two free parameters, of order unity, which 
account respectively for th e shear and bulk vis cosity in 
the Navier- Stokes equation ( Wat kins et al. 19961 ). and for 
the von Neuman-Richtmyer viscosity, which works in high 
Mach number flows. 

A more recent formulation is th e so-called "signal ve- 
locity" viscosity ( Monaghanl Il997 ). which works better 



for small pair separations between particles, being a first- 
order expression of h/r; it reads 



szq 



Pij 



with the signal velocity defined as 



(61) 



(62) 



One or the other of the two formulations above can be 
used in EvoL. 

In general, artificial dissipation should be applied only 
where really necessary. However, it has long been recog- 
nized that artificial viscosity un-physically boosts post- 
shock shear viscosity, suppressing structures in the veloc- 
ity field on scales well above the nominal resolution limit, 
and damping the generation of turbulence by fluid insta- 
bilities. 

To a void the overestimate of the shear viscosity, 
iBalsara (|l995h proposed to multiply the chosen expres- 



sion of liij by a function fij = (/^ + /j)/2; the proposed 
expression for the "limiter" fi for a particle i is 



(63) 



IV, 



^1 + |Vi X Vi\ -^rjCi/hi' 



where rj ~ 10~^ is a factor introduced to avoid numerical 
divergences. It can be seen that / acts as a limiter reducing 
the efficiency of viscous forces in presence of rotational 
flows. 

Another possible cure to the problem is as follows. The 
most commonly used values for the parameters a and P in 
Eq. [571 are = 1 and P = 2. Bulk viscosity is primary pro- 
duced by a, and this sometimes ove r-smoothen the post- 
shock turbulence. To cope with this, iMorris fc Monaghan 



(|l997i) have proposed a modification to Eq. [571 which 
the parameter a depends on the particle and changes with 
time according to: 



dai 
~dt 



Git 



(64) 



where arain = 0.01 is the minimum allowed value, 
r = O.lhi/v^'^^ is an e-folding time scale (with v^'^^ = 
mdi-Xj[v^j^]), and Gi is a source term which can be pa- 
rameterized as Gi = 0.75/i max[0, — |V • Vi\]. This for- 
mulation embodies some different prescriptions (e.g. in 
Rosswog fc Pried I2QQ7I : |Price"2008'). In Eq. [Ml one can 
then put a = ^{c^i + c^j). [Dolag et al. (2005) have shown 
that this scheme captures shocks as well as the original for- 
mulation of SPH, but, in regions away from shocks, the nu- 
merical viscosity is much smaller. In their high resolution 
cluster simulations, this resulted in higher levels of turbu- 
lent gas motions in the ICM, significantly affecting radial 
gas profiles and bulk cluster properties. Interestingly, this 
tends to reduce the differences between SPH and adaptive 
mesh refinement simulations of cluster formation. 
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3.3.2. Artificial thermal conductivity 



As [Price (2008) pointed out, an approach similar to that 
described above should be used to smooth out disconti- 
nuities in all physical variables. In particular, an artificial 
thermal conductivity is necessary to resolve discontinuities 
in thermal energy (even if only a few formulations of SPH 
in literature take this into account). 

The artificial thermal conductivity is represented by 
the term 



Pij 



(65) 



Here is a conductivity coefficient, ^(a^+ap that varies 
with time according to 



~dt 



QU 



(66) 



where r is the same time scale as above, ol^ij^ is usually 0, 
and the source term can be defined as S'f = hi\V'^Ui\/ ^^/ul 
) , in which t he expression for the second deriva- 



(|Pricell2C 

tive is taken from iBrookshaw ( 19861 ) 



V'^Ui ^ 2m^ 



Ui - Uj \ViWij\ 



Pj 



(67) 



which reduces particles noise; v^^^ can be either the same 
quantity defined in Eq. [62] or 



Pij 



(68) 



In passing, we note that adopting the source term sug- 
gested by Price fc M onaghan (2005), S]" = OAhi\\/^Ui\, 
would be less accurate in problems requiring an efficient 
thermal conduction such as the standard shock tube test 
with un-smoothed initial conditions, see Sect. [5) 

The term 11^ must finally be added to Eq. [57J giving 



dui 
~dt 



(69) 



E 



^1 



{vj-Vi)-ViWij{hi) 



As in the case of artificial viscosity, it is worth noting that 
this conductive term is not intended to reproduce a phys- 
ical dissipation; instead, it is a merely numerical artifact 
introduced to smooth out unresolvable discontinuities. 



3.4. Entropy equation 

Instead of looking at the thermal energy, we may follow 
the evolution of a function of the total particle entropy, 
for instance the entropic function 



Ais) 



07' 



(70) 



where 7 is the adiabatic index. This function is expected 
to grow monotonically in presence of shocks, to change 
because of radiative losses or external heating, and to 
remain constant for a n ad iabatic evolution of the gas. 
[Springel fc HernquistI ([200 2) suggest the following SPH 
expression for the equation of entropy conservation: 



dAj 
dt 



7 - 1 

—S{pi,Ui) 



(71) 



2 p] 



where S is the source function describing the energy losses 
or gains due to non adiabatic physics apart from viscos- 
ity. If S << 1, the function A{s) can only grow in time, 
because of the shock viscosity. Eq. [72] can be integrated 
instead of Eq. [57] giving a more stable behaviour in some 
particular cases. Note that internal energy and entropic 
function of a particle can be related via 



Ais) 



(72) 



The entropic function can be used to detect situations of 
shocks, since its value only depends on the strength of the 
viscous dissipation. 

In general, if the energy equation is used (integrated), 
the entropy is not exactly conserved, and viceversa. 
Moreover, if the density is evolved according to the conti- 
nuity Eq. [38] and the thermal energy according to Eq. [SH 
both energy and entropy are conserved; but in this case 
the mass is not be conserved (to ensure this latter, den- 
sity must be computed with Eq.[34j). This problem should 
be cured with the inclusion of the \/h terms in the SPH 
method, as already discussed. 

3.5. X-SPH 

As an optional alternative, EvoL may adopt the smooth- 
ing of the velocity field by replacing the normal equation 
of motion dvi/dt = Vi with 



dvi 
~dt 



^E 



P'ji 



(73) 



( Monaghanl I1992I ) with pji = {pi + pj)/'^^ Vji = Vj — Vi 



and T] constant, < 77 < 1. In this variant of the SPH 
method, known as X-SPH, a particle moves with a veloc- 
ity closer to the average velocity in the surroundings. This 
formulation can be useful in simulations of weakly incom- 
pressible fluids, where it keeps particles in order even in 
absence of viscosity, or to avoid undesired inter-particle 
penetration. 

4. Description of the Code: Integration 

4.1. Leapfrog integrator 

Particles positions and velocities are advanced in time by 
means of a standard leapfrog integrator, in the so-called 
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Kick- Drift- Kick (KDK) version, where the K operator 
evolves velocities and the D operator evolves positions: 



K1/2: 

D: rri+i = Vn 

K1/2: Vn-\-l 



-a{rn)St 



■ V. 



ri+l/2 



St 



1+1/2 



-a{rn-\-i)St. 



A predicted value of physical quantities at the end of each 
time-step is also predicted at the beginning the same step, 
to guarantee synchronization in the calculation of the ac- 
celerations and other quantities. In particular, in the com- 
putation of the viscous acceleration on gas particles the 
synchronized predicted velocity Vn+i = Vn -\- a{rn)St is 
used. 

Moreover, if the individual time-stepping is adopted 
(see below. Sect. I4.2.T]) . non-active particles use predicted 
quantities to give their own contributions to forces and 
interactions. 

4.2. Time-stepping 

A trade-off between accuracy and efficiency can be reached 
by a suitable choice of the time- step. One would like to 
obtain realistic results in a reasonable amount of CPU 
time. To do so, numerical stability must be ensured, to- 
gether with a proper description of all relevant physical 
phenomena, at the same time reducing the computational 
cost keeping the time-step as large as possible. 

To this aim, at the end of each step, for each particle 
the optimal time-step is given by the shortest among the 
following times: 



^tacc,i — Va 




with obvious meaning of symbols; the parameters 77 's are 
of the order of 0.1. The third of these ratios, under par- 
ticular situations, can assume extremely small values, and 
should therefore be used with caution. 

We can compute two more values, the ffist obtained 
from the well known Courant condition and the other con- 
structed to avoid large jumps in thermal energy: 



hi 



Ui 



dui I dt 



where a, (5 and /i are the quantities defined in the artifi- 
cial viscosity parametrization, see Sect. 13.3. If and c is the 
sound speed. 

If the simulation is run with a global time-stepping 
scheme, all particles are advanced adopting the minimum 
of all individual time-steps calculated in this way; syn- 
chronization is clearly guaranteed by definition. 

If desired, a minimum time-step can be imposed. In 
this case, however, limiters must be adopted to avoid nu- 
merical divergences in accelerations and/or internal ener- 
gies. This of course leads to a poorer description of the 
physical evolution of the system; anyway, in some cases a 
threshold may be necessary to avoid very short time-steps 
due to violent high energy phenomena. 

4.2.1. Individual time-stepping 

Cosmological systems usually display a wide range of den- 
sities and temperatures, so that an accurate description of 
very different physical states is necessary within the same 
simulation. Adopting the global time-stepping scheme may 
cause a huge waste of computational time, since gas parti- 
cles in low-density environments and collisionless particles 
generally require much longer time-steps than the typical 
high-density gas particle. 

In EvoL an individual time-stepping algorithm can 
be adopted, which makes a better use of the computa- 
ti onal resources. It is based on a power s-of- two scheme, as 
in Hernquist fc Katz (|l989l ). All individual particles real 



time steps, Sti^true^ are chosen to be a power-of-two sub- 
division of the largest allowed time-step (which is fixed at 
the beginning of the simulation), i.e. 



St,, 



i,true 



St 72"^* 



(74) 



where rii is chosen as the minimum integer for which the 
condition Sti^tme ^ St, is satisfied. Particles can move to 
a smaller time bin (i.e. a longer time step) at the end 
of their own time step, but are allowed to do so only if 
that time bin is currently time-synchronized with their 
own time bin, thus ensuring synchronization at the end of 
every largest time step. 

A particle is then marked as "active" if ti -\- Sti^true = 
tsys^ where the latter is the global system time, updated 

as tsys,new — ^sys^old H~ Stgyg^ where Stgys — TfliTii^Sti) . At 

each time-step, only active particles re-compute their den- 
sity and update their accelerations (and velocities) via 
summation over the tree and neighboring particles. Non- 
active particles keep their acceleration and velocity fixed 
and only update their position using a prediction scheme 
at the beginning of the step. Note that some other evo- 
lutionary features (i.e. cooling or chemical evolution) are 
instead updated at every time-step for all particles. 

Since in presence of a sudden release of very hot ma- 
terial (e.g. during a Supernova explosion) one or a few 
particles may have to drastically and abruptly reduce 
their time-step, the adoption of individual time-stepping 
scheme would clearly lead to unphysical results, since 
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neighbouring non-active particles would not notice any 
change until they become active again. This is indeed a 
key issue, very poorly investigated up to now. To cope 
with this, a special recipe has been added, forcing non- 
active particles to "wake up" if they are being shocked 
by a neighbour active particle, or in general if their time- 
step is too long (say more than 4-8 times bigger) with re- 
spect to the minimum time-step among their neighbours. 
In this way, when an active particle suddenly changes 
its physical state, all of its neighbouring particles are 
soon informed of the change and are able to react prop- 
erly. ISaitoh fc Makinol (|2QQ8l ) recently studied the prob- 
lem, coming to similar conclusions and recommending a 
similar prescription when modelling systems in which sud- 
den releases of energy are expected. 

Note that when a non-active particle is waken up en- 
ergy and momentum are not perfectly conserved, since its 
evolution during the current time-step is not completed. 
Anyway, this error is negligible if compared to the errors 
introduced by ignoring this correction. 

4.3. Periodic boundary conditions 

Periodic boundary c onditions are implemented in EvoL 
using the well known lEwaldl (|l92l[ ) method, i.e. adding to 
the acceleration of each particle an extra-term due to the 
replicas of particles and nodes, expanding in Fourier series 
the potential and the density fluctuation, as descri bed in 
Hernquist et all (|l99lh and in lSpringel et al.l (|2QQll ). If Xi 
is the coordinate at the point of force-evaluation relative 
to a particle or a node j of mass rrij , the additional accel- 
eration that must be added to account for the action of 
the infinite replicas of j is given by 



CLCCpQf(^X^ — 'TTL j 



erfc{uo\x — nL\ 



X ^ 



X — nL 
\x — nL\^ 



(75) 
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27r 
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h • X 



where n and h are integer triplets, L is the periodic 
box size, and uo is an arbitrary number. Convergence is 
achie ved for uo = 2/L, \n\ < 5 and \h\ < 5 (jSpringel et al. 
2QQlh . The correcting terms acCper/TUj are tabulated, and 
trilinear interpolation off the grid is used to compute the 
correct accelerations. It must be pointed out that this in- 
terpolation significantly slows down the tree algorithm, 
almost doubling the CPU time. 

Periodism in the SPH scheme is obtained straightfor- 
wardly finding neighbours across the spatial domain and 
taking the modulus of their distance, in order to obtain 
their nearest image with respect to the active particle. In 
this way, no ghost particles need to be introduced. 



4.4. Parallelization 

EvoL uses the MPI communication protocol to run in 
parallel on multiple CPUs. First of all, particles must be 
subdivided among the CPUs, so that each processor has 
to deal only with a limited subset of the total number of 
bodies. To this aim, at each time-step the spatial domain 
is subdivided using the Orthogonal Recursive Bisection 
(ORB) scheme. In practice, each particle keeps track of 
the time spent to compute its properties during the pre- 
vious time-step in which it has been active, storing it in a 
work- load variable. At the next time-step the spatial do- 
main is subdivided trying to assign to each CPU the same 
total amount of work-load (only considering active par- 
ticles if the individual time-stepping option is selected), 
re-assigning bodies near the borders among the CPUs. To 
this aim, the domain is first cut along the x-axis, at the 
coordinate Xcut that better satisfies the equal work-load 
condition among the two groups of CPUs formed by those 
with geometrical centers Xcentre ^ ^cnt and those with 
Xcentre > Xcut- In practice, cach CPU exchanges particles 
only if its borders are across Xcut- Then, the two groups 
of CPUs are themselves subdivided into two new groups, 
this time cutting along the y-axis (at two different coor- 
dinates, because they are independent from one another 
having already exchanged particles along the x-axis). The 
process is iterated as long as required (depending to the 
number of available CPUS) cutting recursively along the 
three dimensions. It should be noted that this scheme al- 
lows to use only a power-of-two number of CPUs, whereas 
other procedur es (e.g. a Peano-Hilbert domain decompo- 
sition, see e.g. ISpringell [2005) can more efficient exploit 
the available resources. 

It may sometimes happen that a CPU wants to acquire 
more particles than permitted by the dimensions of its 
arrays. In this case, data structures are re-allocated using 
temporary arrays to store extra-data, as described in the 
next Section. 

Apart from the allocation of particles to the CPUs, 
the other task for the parallel architecture is to spread the 
information about particles belonging to different proces- 
sors. This is achieved by means of "ghost" structures, in 
which such properties are stored when necessary. A first 
harvest among CPUs is made at each time step before 
computing SPH densities: at this point, positions, masses 
and other useful physical features of nearby nodes and 
particles, belonging to other processors, are needed. Each 
CPU "talks" only to another CPU per time, first import- 
ing all the data-structures belonging to it in temporary 
arrays, and then saving within its "ghost-tree" structure 
only the data relative to those nodes and particles which 
will be actually used. For example, if a "ghost node" is suf- 
ficiently far away from the active CPU geometrical posi- 
tion so that the gravitational opening criterion is satisfied, 
there will be no need to open it to compute gravitational 
forces, and no other data relative to the particles and sub- 
nodes it contains will be stored in the ghost-tree. 
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Fig. 2. Data structures in EvoL. Void cells are included to allow for new particles to be created and/or imported 
from other CPUs. Ghost structures are reallocated at each time- step. 



A second communication among CPUs is necessary be- 
fore computing accelerations. Here, each particle needs 
to known the exact value of the smoothing and soften- 
ing lengths (which have been updated during the density 
evaluation stage) as well as many other physical values of 
neighbouring (but belonging to different CPUs) particles. 
The communication scheme is exactly the same as before, 
involving ghost structures. It was found that, instead of 
upgrading the ghost-tree built before the density evalu- 
ation, a much faster approach is to completely re-build 
it. 

The dimensions of the ghost arrays are estimated at 
the beginning of the simulation, and every N time-steps 
(usually, N = 100), by means of a "preliminary walk" 
among all CPUs. In intermediate steps, the dimensions of 
these arrays are modified on the basis of the previous step 
evaluation and construction of ghost structures. 

4.5. Data structures 

All data structures are dynamically adapted to the size of 
the problem under exploration. In practice, at the begin- 
ning of a simulation, all arrays are allocated so that their 
dimension is equal to the number of particles per proces- 
sor, plus some empty margin to allow for an increase of 
the number of particles to be handled by the processor 
(see Fig. [2]). 

Anyway, whenever a set of arrays within a processor 
becomes full of data and has no more room for new parti- 
cles, the whole data structure of that processor is re-scaled 
increasing its dimensions, to allow for new particles to be 
included. Of course, also the opposite case (i.e. if the num- 
ber of particles decreases, leaving too many void places) 
is considered. 



To this aim, the data in excess are first stored in tem- 
porary T arrays. Then, the whole "old" data structure 
is copied onto new temporary N arrays, of the new re- 
quired dimensions. Finally, T arrays are copied into the 
N structure; this substitutes the structure, which is 
then deleted. Clearly, this procedure is quite memory- 
consuming. Anyway, it ideally has to be carried out only 
a few times per simulation, at least in standard cases. 
Moreover, it is performed after the "ghost" scheme de- 
scribed above has been completed, and ghost-structures 
have been deallocated, leaving room for other memory- 
expensive routines. 

5. Tests of the code 

An extended set of hydrodynamical tests have been per- 
formed to check the performance of EvoL under different 
demanding conditions: 



Rare faction 
19911 ID); 



and expansion problem ([Einfeldt et al 



Shock tube problem ( SodI 



9781. both in ID and in 3D); 



— Interacting blast waves ([Woodward fc Colellal 11984 
ID); 

— Strong shock collapse (lNohlll987l . 2D); 

— Kelvin-Helmholtz instabili ty problem (2D) : 

— Gresho vorticity problem ( Gresho fc Chaniri99d 2D); 

— Point-like explosion and blast wave (Sedov- Taylor 
problem, 3D); 

— Adiabatic collapse (lEvrardlll988l. 3D): 

— Isothermal collapse ( Boss fc Bodenheimerl 19791 3D); 

— Collision of two politropic spheres (3D); 

— Evolution of a two-component fiuid (3D). 

Except where explicitly pointed out, the tests have 
been run on a single CPU, with the global time-stepping 
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algorithm, and adopting the following parameters for the 
SPH scheme: r] 1.2, amax 2, amin 0.01, p 2a, 
7 = 5/3, and equation of state P = (7 — l)pu. Where 
gravitation is present, the opening angle for the tree code 
has been set = 0.8. The tests in one and two dimensions 
have been run adopting the kernel formulations shown in 
Appendix A. The exact analytical solutions have been ob- 
tained using the software facilities freely available at the 
website http:/ / cococubed. asu. edu/code-pages /codes, shtml. 



been started from un-smoothed initial conditions, leaving 
the code to smear out the initial discontinuities. 

5.2.1. 1-D tests 

The first group of tests have been performed in one di- 
mension. 1,000 equal mass particles have been displaced 
on the X axis, in the domain < x < 1, and changing 
inter-particle spacing, to obtain the following setting: 



5.1. Rarefaction wave 

We begin our analysis with a test designed to check the 
the code in extreme, potentially dangerous low-density re- 
gions, where some iterative Riemann solvers can return 
negative values for pressure /density. In the Einfeldt rar- 
efaction test (Einfe ldt et al. 1991) the two halves of the 
computational domain < x < 1 move in opposite direc- 
tions, and thereby create a region of very low density near 
the initial velocity discontinuity at x = 0.5. The initial 
conditions of this test are 



p = 1.0, P = 0.4, = -2.0, for x < 0.5, 
p = 1.0, P = 0.4, V = 2.0, for x > 0.5. 



(76) 



The results at t = 0.2 for a 1-D, 1000-particle calcula- 
tion are shown in Fig. [3l Clearly, the low-density region is 
well described by the code. 



p = 1.0, u = 1.5, V = 0, for X < 0.5, 
p = 0.25, u = 1.077, V = 0, for x> 0.5. 



(77) 



(m = 6.25 X 10~^). In this way, particles belonging to 
the left region of the tube initially have P = 1, whereas 
particles on the right s ide of the discontinuit y have P = 
0.1795, as in the classic iMonaghan fc Gingoldl (|l983h test. 
We performed four runs: 

— Tl - standard viscosity (with a = 1 and P = 2) and 
no artificial conduction; 

— T2 - standard viscosity plus artificial conduction; 

— T3 - viscosity with variable a plus artificial conduction; 

— T4 - the same as T2, but with the same density dis- 
continuity obtained by using particles with different 
masses on a regular lattice instead of equal mass par- 
ticles shrinking their spacing. To do so, particles in the 
left part of the tube have the mass m = 10~^, whereas 
those in the right part have m = 2.5 x 10~^. 




0.5 1 0.5 1 



X X 

Fig. 3. Left to right, top to bottom: density, internal en- 
ergy, X- velocity, and pressure profiles at t=0.2 for the 
Einfeldt test. All quantities are in code units. 



5.2. Shock Tube 

The classic Riemann problem to check the shock capt uring 
capability is the Sod's shock tube test (|Sodl Il978l ). We 
run many different cases changing numerical parameters, 
to check the response to different settings. All tests have 



1.8 




Fig. 4. Left to right, top to bottom: density, internal en- 
ergy, X- velocity and pressure profiles at t=0.12 for test Tl 
see text for details. Red solid line: exact solution. Alia 
quantities are in code units. 

Figs, m through [3 show the density, internal energy, x- 
velocity and pressure x-profiles at t = 0.12 for the four 
cases (in code units). 

The overall behaviour is good in all cases. We note 
the typical blip in the pressure profile and the wall-like 
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Fig. 5. the same as in Fig. HI but for the test T2 (see text Fig. 8. The same as in Fig. HI but for the test T5 (see text 



for details). 



for details). 
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Fig. 6. The same as in Fig. HI but for the test T3 (see text Fig. 9. The same as in Fig. HI but for the test T6 (see text 



for details). 



for details). 
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overshoot in thermal energy at the contact discontinuity 
in the Tl run. The introduction of the artificial thermal 
conductivity (T2) clearly cures the problem at the expense 
of a slightly shallower thermal energy profile. Introducing 
the variable a-viscosity (T3), the reduced viscosity causes 
post-shock oscillations in the density (and energy), and 
makes the velocity field noisier and turbulent, as expected. 

Little differences can be seen in the T4 run, except for 
a slightly smoother density determination at the central 
contact discontinuity, which also makes the blip problem 
worse. Thus, the best results are obtained with the T2 
configuration. 

5.2.2. 3-D tests 



Fig. 7. The same as in Fig. HI but for the test T4 (see text As pointed out bv ISteinmetz fc Muellerl (|l993[ ). perform- 



for details). 



ing shock tube calculations with 1-D codes, with particles 
initially displaced on a regular grid, makes the effects of 
numerical diffusivity essentially negligible, but in realis- 
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tic cosmological 3-D siraulations the situation is clearly 
more intrigued. To investigate this issue, we switch to a 
3-D description of the shock tube problem. The initial 
conditions are set using a relaxed glass distribution of 
~ 60, 000 particles, in a box of dimensions [0, 1] along 
the X axis and (periodic) [0,0.1] along the y and z axis. 
Particles are assigned a mass rriright = 1-71 x 10~^ and 
mieft = 2.5 X rriright on the two sides on the discontinuity 
at X = 0.5, with all other parameters as in the 1-D case. A 
first run (T5) was performed adopting the Tl configura- 
tion (standard viscosity, no artificial conduction), while a 
second one (T6) was performed switching on the artificial 
conduction term. 

Figs. [8] and [9] show the state of the system at t = 0.12 
in the two cases (note that all particles and not mean val- 
ues are plotted). As in the 1-D case, the overall agreement 
with the theoretical expectation is quite good. However, 
the inclusion of the conductive term on one hand yields 
sharper profiles and reduces the pressure blip, on the other 
hand gives a smoother description of contact discontinu- 
ities (see the density profile at x = 0.55). We argue that 
part of the worse results may be due to the use of particles 
with different mass as a similar result has been found with 
the corresponding 1-D test (see T4 above). A limit on the 
conduction efficiency could give better results, but and ad 
hoc adjustment should be applied to each case. 

5.3. Interacting blast waves 

Woodward fc Colellal dlOsl proposed a severe test m 1-D, 



i.e. two strong blast waves develop out of a gas of density 
p = 1 initially at rest in the domain < x < 1, with 
pressure P = 1000 for x < 0.1, P = 100 for x > 0.9 and 
P = 0.01 in between (the adiabatic index is 7 = 1.4). The 
boundary conditions are reflective on both sides; this can 
be achieved introducing ghost particles created on-the-fly. 
Evolving in time, the two blast waves give rise to multiple 
interactions of strong shocks, rarefaction and reflections. 

We performed two runs: one at high resolution with 
1000 equally spaced particles, and one at low resolution 
with only 200 particles. In both cases, the code included 
the artificial conduction, and viscosity with constant a = 
1. 

Fig. [To] shows the density x-profiles for the high reso- 
lution case, at t = 0.010, 0.016, 0.026, 0.028, 0.030, 0.032, 
0.034, 0.038. The results can be compared with those o f 
Woodward Colellal ([TqsI . LSteinmetz Muellej (|l993h 
and ISpringell (|2009l ). The complex evolution of the shock 
and rarefaction waves is reproduced with good accuracy 
and with sharp profiles, in nice agreement with the the- 
oretical expectations. The very narrow density peak at 
t = 0.028, which should have p ^ 30, is well reproduced 
and only slightly underestimated. 

Fig. [m shows the comparison between the density pro- 
files of the high and low resolution cases, at t = 0.038. The 
low-resolution one shows a substantial smoothing of the 
contact discontinuities; nevertheless, the overall behaviour 



is still reproduced, shock fronts are located at the correct 
positions, and the gross features of the density profile are 
roughly reproduced. 

5.4. Kelvin-Helmoltz instability 

The Kelvin-Helmoltz (K H) instabil ity problem is a very 
de manding test. Recen tly, IPricd ( 2008.) replied to the claim 
bv lAgertz et al ] (|2007h that the SPH algorithm is not able 
to correctly model contact discontinuities like the Kelvin- 
Helmoltz (KH) problem. He showed how the inclusion of 
the conductive term is sufficient to correctly reproduce 
the instability. We try to reproduce his results running a 
similar test, using the 2-D version of EvoL. 

We set up a regular lattice of particles in the periodic 
domain 0<x<l,0<y<l, using 256^ particles . The 
initial conditions are similar to those adopted by IPrice 
(2008). The masses of particles are assigned in such a way 
that p = 2 fo r \y — 0.5 1 < 0.25, and p = I elsewhere 
(note that in IPrice 20081 this density difference was in- 
stead obtained changing the disposition of equal mass par- 
ticles). The regions are brought to pressure equilibrium 
with P = 2.5. In this case, although the initial gradi- 
ents in density and thermal energy are not smoothed, the 
thermal energy is assigned to the particles after the first 
density calculation, to ensure a continuous initial pressure 
field. 




Fig. 11. Comparison between high-resolution (solid line) 
and low-resolution (dots) density profiles at t = 0.038 for 
the Woodward test (see text for details). The quantities 
are in code units. 



The initial velocity of the fluid is set as follows: 



[0.5 if |^-0.5| < 0.25, 
I -0.5 if |^-0.5| > 0.25, 



(78) 



and 
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Fig. 10. Left to right, top to bottom: density profiles at t = 0.010, t = 0.016, t = 0.026, t = 0.028, t = 0.030, t = 0.032, 
t = 0.034, t = 0.038 for the Woodward high-resolution test (see text for details). All quantities are in code units. 
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where we use a = 0.05/a/2, and wq = 0.1 or wq = 0.25 for 
two different runs (Kl and K2, respectively). In this way, 
a small perturbation in the velocity field is forced near the 
contact discontinuity, triggering the instability. 

Figs. [12] and [13] show the temporal evolution of the 
instability in the two cases. Clearly, our code is able to 
reproduce the expected trend of whirpools and mixing 
between the fluids; the evolution is faster for a stronger 
initial perturbation. For comparison, we also run a test 
with the K2 initial conditions, but switching off artificial 
conduction; Fig. [14] shows the comparison between the fi- 
nal density fields (at t = 1.79), and between the initial 
pressure fields. In the latter, the blip that inhibits mixing 
is clearly visible at the discontinuity layer: the resulting 
spurious surface tension keeps the two fluids separated, a 
problem efficiently cured by the artificial conduction term. 

5.5. Strong shock collapse 

The strong shock collapse (jNohlll987 ) is another very de- 
manding test, and is generally considered a serious chal- 
lenge for AMR numerical codes. We run the test in 2-D, 
displacing ^ 125,000 particles on a regular grid and re- 
tailing those with r < 2 to model a gas disk, with ini- 
tially uniform density p = 1, and negligible internal en- 
ergy; we then add a uniform radial velocity towards the 
origin, Vin = — 1. A spherical shock wave with formally 
infinite Mach number develops from the centre and trav- 
els outwards, leaving a constant density region inside the 
shock (in 2-D, pint = 16). The test was run on 4 parallel 
CPUs. 



In Fig. [151 the density in the x > 0, y > region is 
shown at t = 1.2, while in Fig. [16] the density-radius rela- 
tion at the same time is plotted for all particles. There is 
good agreement with the theoretical expectation, although 
a strong scatter is present in particles density within the 
shocked region, where moreover the density is slightly un- 
derestimated (anywa y, this is a co mmon feature in this 
kind of test: see e.s;. ISDringell[2009h . 

Apart from these minor flaws, the code is able to han- 
dle the very strong shock without any particular difficulty. 

5.6. Gresho vortex problem 

The test for the conservation of angular momentum 
and vorti city, the so-call ed Gresho triangular vortex 
problem (ICresho fc Chan [l^90) . in its 2-D version 
(|Liska fc Wendrofjl2003l ). follows the evolution of a vor- 
tex initially described by an azimuthal velocity profile 



5r if0<r<0.2, 
2-5r if0.2<^i<0.4, 
if r > 0.4, 



(80) 



in a gas of constant density p = 1. The centrifugal force is 
balanced by the pressure gradient given by an appropriate 
pressure profile. 



P(r) = 
5 + 12.5r2 
9 + 12.5r2 
3 + 41n(2) 



(81) 



if < r < 0.2, 
20r + 41n(5r) if 0.2 <u< 0.4, 
if r > 0.4, 



so that the vorte x should remai n independent of time. 

As noted bv ISpringell ( 20091 ). this test seems to be a 
serious obstacle for SPH codes. Due to shear forces, the 
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Fig. 12. Kelvin-Helmoltz instability according to test Kl. Left to right, top to bottom: density field at t = 0.768, 
t = 1.158, t = 1.549 and t = 1.940. All quantities are in code units. 
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Fig. 13. Kelvin-Helmoltz instability according to test K2. Left to right, top to bottom: density field at t = 0.377, 
t = 0.767, t = 1.158 and t = 1.549. All quantities are in code units. 



angular momentum tends to be transferred from the in- the rotational motion to spread within the r > 0.4 re- 
ner to the outer regions of the vortex, eventually causing gions; these finally dissipate their energy because of the 
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Fig. 14. Kelvin- Helmoltz instability according to test K2, with (left) and without (right) artificial thermal conductivity. 
Top panels: initial pressure field. Bottom panels: density fields at t = 1.79. All quantities are in code units. 



interactions with the counter-rotating regions of the peri- 
odic replicas. With the GadgetS code, the vortex did not 
survive up to t = 3. AMR codes generally give a better 
description, although under some conditions (for example, 
the addition of a global bulk velocity, which should leave 
the system unchanged due to its galileian invariance) they 
also show some fiaws. 

We tried to run the test starting from two differ- 
ent initial s etting s of particles, since as pointed out by 
Monaghanl (|2QQ6l ) the initial particle setting may have 
some infiuence on the subsequent evolution of the sys- 
tem. Thus, in the first case we set up a regular carte- 
sian lattice of ~ 10, 000 particles (in the periodic domain 
—0.5 < X < 0.5, —0.5 < y < 0.5), and the second one put 
the same number of particles on concentric rings, allowing 
for a more ordered description of the rotational features 
of the gas (see Fig. [TTj) . 

The results are shown in Fig. [181 We find a resid- 
ual vorticity still present at t ~ 3, although strongly de- 
graded and reduced in magnitude. As expected, a substan- 
tial amount of the angular momentum has been passed to 
the outer regions of the volume under consideration. This 
is confirmed by comparing the initial to the final pressure 
fields (Fig. [19]): at t = 2.7 the pressure exterted by the out- 
skirts is higher, due to the increased temperature of the 
regions confining with the periodic replicas. This leads to 
a compression of the internal region, which in turn must 
increase its pressure as well. 

However, we could not find in literature other direct 
comparisons for this test with an SPH code. It must be 



pointed out that the resolution in these tests is quite low, 
and a larger number of particles may help improving the 
results. 



5. 7. Point explosion 

The presence of the artificial conduction becomes crucial 
in the point explosion experiment, a classic test to check 
the ability of a code to follow the evolution of strong 
shocks in three dimensions. In this test, a spherically 
symmetric Sedov- Taylor blast-wave develops and must be 
modelled consistently with the known analytical solution. 

We run the test in 3-D. Our initial setting of the SPH 
particles is a uniform grid, in a box of size [—0.5, 0.5] along 
each dimension, with uniform density p = 1 and negligi- 
ble initial internal energy. The grid is made of 31^ par- 
ticles. At t = an additional amount £^ = 1 of thermal 
energy is given to the central particle; the system is sub- 
sequently let free to evolve self-consistently (without con- 
sidering self-gravity). It must be noted that, again, the 
central additional thermal energy is not smoothed, so to 
obtain extreme conditions. 

The following test runs are made (using the 
Lagrangian SPH scheme): 

— SI - constant viscosity with a = 1 viscosity, no artifi- 
cial conduction; 

— S2 - variable viscosity with a plus artificial conduction; 

— S3 - as S2, but without the \/h correcting terms; 

— S4 - as S2, with the individual time-stepping scheme; 
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S5 - as S2, with higher resolution (51^ particles). 




Fig. 15. The Noh's problem: density field in the first quad- 
rant at t = 1.2. All quantities are in code units. See text 
for details. 
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Fig. 16. The Noh's problem: density vs. radius relation for 
all particles at t = 1.2. Solid line: theoretical expectation. 
All quantities are in code units. See text for details. 

Fig. [2Q1 shows the xy positions of particles belonging to 
the z = plane at t = 0.09 for the SI test, and the corre- 
sponding radial density profile, compared to the analytical 
expectation. A strong noise (which develops soon after the 
initial explosion) is evident in the particle positions, due to 
the huge unresolved discontinuity in thermal energy. The 
description of the shock evolution is quite poor, although 
its gross features are still captured. 

Fig. [2T] presents the same data for the S2 run, in which 
both the viscosity with variable a and artificial conduction 
are switched on. The increased precision in the descrip- 
tion of the problem is immediately evident: the conductive 
term acts smoothing out the thermal energy budget of the 
central particle, strongly reducing the noise and allowing 
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Fig. 17. Positions of all particles at t = 0.02 (top) and t = 
2.7 (bottom) for the two different initial configurations for 
the Gresho test: grid (left) and rings (right). All quantities 
are in code units. See text for details. 




Fig. 18. Initial (black dots) and final (red small points) 
azimuthal velocity profiles for the two Gresho tests: grid 
(top) and rings (bottom). All quantities are in code units. 
See the text for details. 
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Fig. 19. Initial {t = 0.02, left) and final {t = 2.7, right) 
pressure fields in the grid Gresho test. Note the increased 
pressure in the outskirts (white) and in the central regions 
of the vortex. All quantities are in code units. See the text 
for details. 
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Fig. 20. Top: spatial distribution of particles in a slice 
taken at z=0 and t=0.09 for the Sedov blast wave of 
test SI. Bottom: Radial density profile for the same test. 
Points: SPH particles; solid line: analytical prediction. All 
quantities are in code units. 



a symmetric and ordered dynamical evolution. The radial 
density profile (in which single particles, and not averages, 
are plotted) shows that the shock is well described. The 
density peak is a factor of ~ 2 lower than the analytical 
expectation, due to the intrinsic smoothing nature of the 
SPH technique. 

Another test (S3) is carried dropping the \/h correct- 
ing terms (see Sect. [3j). While the differences in particle 
positions and density profiles are almost negligible, the 
comparison of the total energy conservations E{t)/ E{to) 
in the two cases (Fig. [22]) shows that the correction helps 
in containing losses due to wrong determinations of the 
gradient. 

A fourth test (S4) is also run to check the behaviour 
of the individual time-stepping algorithm; the results are 
virtually indiscernible from the those of the previous case 
(not shown here). This is noteworthy as all particles must 
be "waken up" (using the method described in Sect. 14. 2. T]) . 
being still and cold at the beginning of the simulation. 

Finally, a test (S5) with increased numerical resolution 
is performed (Fig. [23]). As expected, in this case, while the 
particles positions again accurately describe the evolution 



Fig. 21. The same as in Fig. [20l but for the S2 test. 




Fig. 22. Energy conservations for Sedov test with artificial 
conductivity and variable viscosity. Points: with Vh terms; 
crosses: without \/h terms. All quantities are in code units. 



of the system, the shock boundary is sharper, and closer 
to the peak value of the analytical solution. 



5.8. Self-gravitating collapse 

The collapse of a gaseou s self-gravitati ng sphere is another 
standard SPH 3-D test (lEvrardlllQsif ). The combined ac- 
tion of hydrodynamics and self-gravity leads the system, 
a sphere of gas initially far from equilibrium, with negli- 
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Fig. 23. The same as in Fig. [20l but for the test S5. 

gible internal energy and density profile p cx r~^, to a col- 
lapse in which most of its kinetic energy is converted into 
heat. An outward moving shock develops, a slow expan- 
sion follows and, at late times, a core-halo structure de- 
velops with nearly isothermal inner regions and the outer 
regions cooling adiabatically. To set the initial conditions, 
10,000 particles of mass 1.25 x lO^^M© are placed in a 
sphere of radius 5 x 10~^ Mpc, in such a way that their 
initial density radial profile scales as 1/r. Another sphere 
with 40,000 particles and consequently smaller particles' 
masses (higher resolution) is also prepared. Two cases are 
examined: namely the collisionless collapse and the adia- 
batic collapse. 

5.8.1. Collisionless collapse 

To check the performance of the adaptive softening length 
algorithm, tests are run switching off hydrodynamical in- 
teractions to mimic the collisionless collapse. A first case 
drops the hydrodynamical interactions, includes the adap- 
tive softening algorithm but neglects the correcting terms 
presented in Section 12. li It is meant to provide the refer- 
ence case. The second one includes also the effects of these 
correcting terms. 

Fig. [Ml shows, in the upper panels, the spatial distribu- 
tion of particles projected onto the xy plane after during 





2000 4000 6000 8000 1000012000 



2000 4000 6000 8000 1000012000 



Fig. 24. Collisionless collapse with adaptive softening 
lengths. Top: spatial distribution of particles projected 
onto the xy plane after 12,000 steps; middle: softening 
lengths (code units) as a function of the radial coordi- 
nate; bottom: temporal conservation of the total energies 
of the systems E{t)/ E{to)^ for the un-corrected (left) and 
corrected (right) tests. 



the collapse of the sphere, for the un-corrected (left) and 
corrected (right) cases: no evident difference is present. 
The same holds for (middle panels) the value of the soft- 
ening lengths of particles as a function of the radial co- 
ordinate. However, the temporal conservation of the total 
energies of the systems E{t)/E{to) in the two cases are 
quite different (bottom panels of Fig. [24j). If the correct- 
ing terms are not included a secular error in the energy 
conservation develops, reaching 2%, while the error is ~ 10 
times smaller (albeit not completely eliminated) with the 
introduction of the extra-terms. 

5.8.2. Adiabatic collapse 

The standard hydrodynamical tests ( EvrardI 1988[ ) are cal- 
culated with the following settings: 

— El - constant softening lengths (e^ = O.li^x (m^/M)^-^, 
with R and M total radius and mass of the sphere 
at the beginning of the simulation); Lagrangian SPH 
scheme, with constant viscosity (a = 1) and no artifi- 
cial conduction; 

— E2 - as El, except for using adaptive softening lengths; 

— E3 - adaptive softening lengths, viscosity with variable 
a, and artificial conduction; 
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- E4 - as E3, higher resolution (40,000 particles). 



The global energy trends are shown in Fig. [25l while 
Fig. [26] shows a magnification of the energy conservation 
E{t)/E{to) for the El and E2 runs. The energy conserva- 
tion is always good, reaching a maximum error of ~ 0.4% 
at the moment of maximum compression. For comparison, 
in th e same test the finite differences method by Evrard 
( 19881) ensures conservation at :^ 1%, while iHarfst et al.l 
( 20061 ) finds a maximum error of ~ 3% with a standard 
SPH implementation. Almost no difference can be noted 
between the variable and the constant e runs, despite a 
variation of more than two orders of magnitude in the 
softening length values in the E2 run (see Fig. [27j). 
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Fig. 25. Energy trends in the cold collapse test. Light 
lines: red: test El; black: test E2; blue: test E3. Heavy dot- 
ted blue line: test E4 (see text for details). Top to bottom: 
thermal, kinetic, total and potential energies. All quanti- 
ties are in code units. 
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Fig. 26. Magnification of the energy conservations 
E{t)/E{to) in the cold collapse tests. Black: testEl; red: 
testE2. See text for details. 
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Fig. 27. Softening lengths in the cold collapse test T2 at 
t = 0.8. The solid line is the value of the constant softening 
length in test Tl. All quantities are in code units. 



The density and internal energy profiles at t = 0.8, for 
the El to E3 runs, are shown in Fig. [28l the comparison 
between the E3 and the hi-res E4 runs is shown in Fig. 
l29l (the shock is moving outward leaving an inner, heated 
core). The discontinuity in density is smoothed out by the 
kernel smoothing, which in the current implementation re- 
quires a high number of neighbors (~ 60), but the effects of 
the inclusion of the artificial conductive term are evident. 
We didn't find any particular problem with the inclusion 
of the conductive term in a self-gravitating system. 

We also run a X-SPH test, adopting r] = 0.25. While 
the energy trends are essentially identical to the ones 
the previous tests (not shown), the energy conservation 
is slightly worse , reaching an error of ~ 2.5%. Indeed, 
Monaghanl ( 2002h pointed out that a more complex treat- 
ment of the equations of motion should be adopted to 
achieve perfect energy conservation. 



5.9. Isothermal collapse 

Introduce d bv IBoss fc Bodenheimer J 1979h and later re- 
visited bv llBurkert Bodenheimei] ('l993V the test on the 
isothermal collapse couples hydrodynamics and gravity, 
also adding rotation. It develops in a complex game of 
collapse and fragmentatio n. We run th e test in the ver- 
sion proposed by B ate fc BurkertI (|l997l ). 

A homogeneous sphere of cold gas with initial radius 
R = 5x 10^^ cm and mass M = IMq is modelled retailing 
~ 78, 000 particles out of a uniform lattice. The sphere is 
put in solid body rotation around the z axis, with angular 
velocity Q = 7.2 x 10""*^^ rad s~^, and the otherwise flat 
density field is perturbed by 10% in such a way that 



p((/)) =po[l + 0.1co5(2(/))]. 
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Fig. 28. Density (top) and internal energy (bottom) pro- 
files at t = 0.8 in the cold collapse El (left), E2 (center) 
and E3 (right) tests. All particles are plotted. 
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Fig. 29. Comparison between cold collapse E3 (low-res, 
left) and E4 (hi-res, right) test: density (up) and internal 
energy (bottom) profiles at t = 0.8. 



where (j) is the azimuthal angle about the rotation axis 
and po = 3.82 x 10~^^ g cm~^ (we achieve the density 
perturbation by varying the masses of particles instead 
of their positions). The gas is described by an isothermal 
equation of state, P = c^p, with initial sound speed Cg = 
1.66 X 10^ cm s~^. It is assumed to be optically thin so 
that no mechanical heating can occur, and the adiabatic 
index is 7 = 1.4. 

The test is very challenging and has been used to check 
the behaviour of bo th AMR and SPH codes. In particu- 
lar, iBate fc BurkertI (|l997l) used it to study the effects of 
wrong calibrations of the smoothing and softening lengths 
in SPH (see Sect. 13. 2p . We run the test in four different 
cases: 

— Bl - Lagrangian SPH with constant softening lengths 
(e 5.26 X 10^^ cm, an intentionally large value); 



— B2 - Lagrangian SPH with adaptive softening lengths 
(no hmits imposed on e and h); 

— B3 - as B2, but imposing a minimum smoothing and 
softening length hmin = ^mAn = 10^"^ cm, to avoid 
artificial clumping (as in Bate fc Burkert|[l997[ ): 

— B4 - as B2, with the X-SPH method (adopting r] = 
0.25). 

Figs. [30] to [33] plot the temporal evolution of the den- 
sity field in the z = plane; shown are snapshots at 
t = 1.0, 1.15, 1.23, and 1.26, in units of free fall time, 
tff = 1.0774 X 10^2 s. Fig.[3lplots the final density field 
at t = 1.29 for the four tests. These plot s can be compar ed 
to those bv lBate fc BurkertI (|l997h and lSpringel (|2005h . 

The large value of e assumed in the Bl test clearly 
demonstrates how a wrong choice of the softening length 
can result in catastrophic errors. The evolution of the sys- 
tem is clearly differen t from the grid solution used by 
iBate fc BurkertI (|l997h as a proxy of the exact solution of 



the problem. Because of the excessive softening, clumps 
may be stable against collapse and a rotating structure 
forms, resembling a barred spiral galaxy. More realistic 
choices for e may give better results, but the opposite 
problem (artificial clumping of unresolved structures) may 
arise in case of too small a softening length. Indeed, the 
case B2, in which the adaptive algorithm is adopted (note 
that e = /i), gives much better results, but disorderly 
clumped sub-structures may form. 

Imposing a minimum e (=/i, case B3) cures the prob- 
lem (note however that the evolution of the system seems 
to be somewhat faster than the grid solution). Anyway, 
this prescription cannot be generalised. A more physically 
sounded solution may be the introduction of a pressure 
floor in unresolved regions, i.e. where the Jeans mass of 
the gas is not resolved by the SPH technique. This is- 
sue will be discussed in the forthcoming companion paper 
(Merlin et al. 2009, in preparation). 

Finally, the X-SPH test B4 gives results that are es- 
sentially undiscernible from those of case B2, apart from 
a slightly less disordered density field in the regions sur- 
rounding the central collapsed structure. 

Fig. [35] plots the temporal evolution of the maximum 
value of the density. Thi s plot can be compa red with 
Fig. 5 a in lBate fc BurkertI (|l997h and Fig. 12 in lSpringell 



(|2005f ). The results for all tests are very similar and agree 
quite well with their SPH tests of comparable resolution 
(with the plateau at p 2 x 10~^^ g cm~^ reached at 
t = 1.14), except for the case Bl which clearly evolves 
faster. 

For reference. Fig. [36] shows the initial and final values 
of both e and h (which are equal in the cases B2, B3 and 
B4) as a function of p in the four runs. 

5.9.1. Parallel run 

Finally, we re-computed the case B2 test using 16 par- 
allel CPUs. Fig. [37] shows the comparison between the 
final density field for this run and for the original one 
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Fig. 30. Evolution of the density field in the central region Fig. 33. The same as in Fig. [30l but for test B4 
for the test Bl. Plotted is the region within 1.54 x 10"^^ 
cm from the origin, on the xy plane. Left to right, top to 
bottom: t = 1.0, 1.15, 1.23, 1.26, in units of free fall time, 
tff = 1.0774 X 10^2 
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Fig. 31. The same as in Fig. [30l but for test B2. 



J 

— \ — ^ — \ — ^ — \ — ^ 




■ J 





Fig. 34. Final density field at t = 1.29 (in units of free fall 
time) in the central region. Left to right, top to bottom: 
tests Bl, B2, B3, B4. See text for details. 
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Fig. 32. The same as in Fig. [30l but for test B3. 
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Fig. 35. Maximum density in the isothermal collapse 
tests. Bl: red dashed line; B2: green solid line; B3: blue 
thick dots: B4: black dot-dashed line. See text for details. 
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Fig. 36. Initial (blue) and final (red) smoothing length h 
of all particles in the four isothermal collapse tests (left 
to right, top to bottom: Bl, B2, B3, B4). The softening 
length e is equal to h in the B2, B3 and B4 runs; the black 
dots in the top left panel show its constant value in the 
test Bl. The abrupt truncation in the bottom left panel 
is due to the imposed hmin (= ^min) in the test B3. 



with single CPU. Although the agreement is good, and 
no sign of a multi-processor subdivision of the domain is 
present, some differences in local morphological features 
are present; these are more evident looking at Fig. [38l 
where the comparison is made for the smoothing lengths 
values. While the gross features and the minimum and 
maximum values are very similar, differences in the local 
clumping of particles are clearly present. 

We ascribe this flaw to the parallel Oct-Tree scheme. 
The global tree constructed by > 1 CPUs is slightly 
different from the one built by a single processor on the 
same spatial domain: the regions in which two or more 
different CPUs interact are described by a different cell 
architecture depending on N (note that the reason for 
this is that the cells have to be always cubic; indeed, the 
problem should not occur if a binary tree, in which cells 
have adaptive shapes, was used). This causes a slightly 
different approximation of the gravitational force, due to 
the opening cell criterion, giving in turn a slightly different 
dynamical evolution. 

To prove the validity of this speculation, we run again 
the same tests reducing the opening angle ^ to 0.1 instead 
of the standard 0.8. In this way, the approximation given 
by the adoption of the tree structure should become of 
negligible importance (at the expenses of a much longer 
computation time). Indeed, looking at Figs. [39] and [40l we 
can notice an improvement in that the density fields are 
more similar, although still not exactly identical. An ideal 
run with ^ = (i.e., exact particle-particle interaction) 
should give almost identical distributions of particles. 

We must point out that this flaw may be of non- 
negligible importance when dealing with phenomena that 
directly depend on the local density field, such as the star 



formation and feedback processes. The same simulation 
run on a different number of CPUs may give slightly dif- 
ferent results because of this issue. 
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Fig. 37. Comparison between the final density fields for a 
single CPU (left) and a 16 parallel CPUs (right) runs of 
the B2 test. See the text for details. 
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Fig. 38. Comparison between the final /i-radius relation, 
for a single CPU (left) and a 16 parallel CPUs (right) runs 
of the B2 test. See the text for details. 
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Fig. 39. Comparison between the final density fields for a 
single CPU (left) and a 16 parallel CPUs (right) runs of 
the B2 test, in a run with tree opening angle 9 = 0.1. See 
the text for details. 



5.10. Collision of two gas spheres 

The collision between two gaseous Plummer spheres is 
very crucial test for the energy conservation. The tests 
are run in 3-D, with adaptive softening lengths, variable 
a viscosity, and thermal conduction. 
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Fig. 40. Comparison between the final /i-radius relation, 
for a single CPU (left) and a 16 parallel CPUs (right) runs 
of the B2 test, in a run with tree opening angle = 0.1. 
See the text for details. 




Each of the two spheres (of unit mass, unit radius, 
and negligible initial internal energy) are set up randomly 
selecting 10,000 particles in a way that their initial density 
radial profile resembles the Plummer profile, p ex (1 + 
^^-5/2 xheir centres are on the x axis. 

In a first run (CI) we follow a head-on collision, as- 
signing a relative velocity of 1.5 in the x direction to 
the spheres. In a second run (C2) we added a shear 
component, giving relative velocities \/^Vx\ = 0.15 and 
\/\vy\ = 0.075. 



According to lHernquistI ( 19931 ) the violent head-on col- 
lision of two polytropic stars does not conserve the en- 
ergy within ^ 1Q%- Mo dern formulations of SPH (see e.g. 
Rosswog fc Pricell 2007j) can handle the problem with much 



more accuracy, reducing the error below 0.1%. 




Fig. 42. Top: energy trends in the collision test CI; black 
dots: total energy, blue dot-dashed line: potential, red 
dashed line: kinetic, green solid line: thermal. Bottom: 
magnification of the total energy conservation E{t)/ E{to). 



Figs. [4T1 and [43l plot the density and velocity fields at 
different times, on a slice taken at z = 0. In our tests, the 
energy conservation is good (see Figs. l42l and [44j) . with a 



Fig. 44. Top: energy trends in the collision test C2; black 
dots: total energy, blue dot-dashed line: potential, red 
dashed line: kinetic, green solid line: thermal. Bottom: 
magnification of the total energy conservation E{t)/E{to). 

maximum error of ~ 0.2% in the head-on collision, and 
~ 1% error when the shear component is added. 

5.11. Two-components evolution 

As a final test, we studied the evolution of a two- 
components fluid in order to investigate the behaviour of 
the adaptive smoothing and softening length algorithm in 
presence of more than one material. We set up two spheres 
of the same radius, R = 10^ pc, and mass, M = 4xlO^M0, 
using for each sphere the particle distribution adopted adi- 
abatic collapse test by Evrard. One of the two systems was 
considered as made of collisionless bodies, with negligible 
pressure, thus mimicking the behaviour of a Dark Matter 
halo; the other one was instead considered as made of gas, 
with an initial temperature of 1000 K. The system was 
the let free to evolve under the action of self-gravity and 
hydrodynamical interactions, in four different runs with 
the following settings: 

— TCI - constant softening lengths, = O.IR x 
(m^/M)^-^; null initial velocities 

— TC2 - adaptive softening lengths, null initial velocities 

— TC3 - as TC2, with an initial solid body rotation Q = 
2.5 X 10~^^ rad s~^, for both components 

— TC4 - as TC3, but with counter-rotating initial tan- 
gential velocities (same magnitude). 

While the collisionless component soon collapses under 
self-gravity, later reaching relaxation, the gas undergoes a 
phase of expansion because of its internal pressure, but in 
the central region, where it also interacts with collapsing 
Dark Matter, reaching very high densities. 

The free expansion of a system of particles is by itself 
a demanding test, and a loss of total energy is a common 
flaw for Lagrangian codes. Furthermore, the rotation im- 
posed in the last two runs makes the problem more com- 
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Fig. 41. Evolution of the density field in the z = plane for the collision test CI. Left to right, top to bottom: t = 0.05, 
0.15, 0.2, 0.25, 0.3, 0.5. Superposed is the velocity field. 




Fig. 43. Evolution of the density field in the z = plane for the collision test C2. Left to right, top to bottom: t = 0.03, 
0.38, 0.77, 1.16, 1.55, 1.94, 2.72, 3.5, 4.29. Superposed is the velocity field. 



plicated. In our tests, a further complication was given 
by the dynamical interaction within the central regions 
of the system, where the hydrodynamical pressure of the 
gas struggled against self-gravity and the (stronger) grav- 
itational action of the collisionless component. Moreover, 
when the adaptive e scheme was used, gas particles had 
to cope with a large variation of their interaction sphere, 
due to the presence of Dark Matter particles in the cen- 
tral regions, and to their absence in the outskirts, where 
in addition gas particles are far away from one another. 

Fig. [35] plots the energy trends for all the four tests. 
The conservation is very good in all cases, with a max- 
imum error well below 1%; quite surprisingly, the worst 
performance is given by the constant e scheme, while the 
corresponding test with adaptive softening performs bet- 
ter staying below a ~ 0.5% error. 

Fig. [46] shows the evolution of the Dark Matter com- 
ponent in the test TC3, projected on the xy and on the 
xz planes (the case TC4 run gives almost undistinguish- 



able plots). Clear is the early formation of a disk due to 
the initial rotational velocity, and a subsequent expan- 
sion which leaves a dense core with a loosely disky shape. 
On the other hand, the gas component expands with al- 
most radial motion (Fig. [47]) . but in the central regions 
the interaction with the collisionless component plays an 
important role. Indeed, looking at Fig. [48] where a mag- 
nification of the velocity field in the central region at late 
times is shown for the tests TC3 and TC4, the difference 
in the two cases is clear: in the counter-rotation run the 
gas in the innermost zone has inverted its sense of motion 
and co-rotates with the Dark Matter, creating a complex 
pattern of velocities. 

The conservation of angular momenta in the tests TC3 
and TC4 is shown in Fig. [49] We find a maximum error of 
~ 6% in the relative error for the total momentum in the 
counter-rotation test; in this case the absolute violation is 
anyway small because of the very small total initial mo- 
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Fig. 45. Two components test, energy conservation 
E{t)/E{to) for the four runs: dashed red hue, TCI; blue 
dot-dashed line, TC2; green dots, TC3; black solid line, 
TC4. 
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Fig. 46. Two components test, case TC3. Dark Matter 
particles projected onto the xy (top panels) and xz (bot- 
tom panels) planes: left to right, top to bottom, t = 0.5, 
1.0, 1.5, 2.0, 2.5 Gyr. 



mentum. In all the other cases the conservation is always 
very satisfactory. 
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Fig. 47. Two components test, case TC3. Gas density and 
velocity fields in a slice taken at z = 0: left to right, top 
to bottom, t = 0.5, 1.0, 1.5, 2.0, 2.5 Gyr. 
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Fig. 48. Two components test: comparison between the 
central velocity fields in cases TC3 (left) and TC4(right), 
in a slice taken at z = 0, at t = 2.5 Gyr. 



6. Conclusions 

We have presented the basic features of EvoL, the new re- 
lease of the Padova parallel N-body code for cosmological 
simulations of galaxy formation and evolution. In this first 
paper, the standard gravitational and hydrodynamical al- 
gorithms have been extensively reviewed and discussed, as 
well as the parallel architecture and the data structures of 
the code. EvoL includes some interesting options such as 
adaptive softening lengths with self-consistent extra-terms 
to avoid large errors in energy conservation, V/i terms in 
the SPH formalism, variable a viscosity formulation, and 
artificial thermal conduction terms to smooth out pressure 
at contact discontinuities. 

We have also performed and presented an extended 
series of standard hydrodynamical tests to check the abil- 
ity of the code to handle potentially difficult situations. 
The results are encouraging. Almost all tests have given 
results in nice agreement with theoretical expectations 
and previous calculations with similar codes, sometimes 
even showing better performance. In particular, we showed 
how the inclusio n of an artific ial thermal conduction term 
as suggested by 'Pricd (l2008l) significantly improves the 
modeling of demanding problems such as the Kelvin- 
Helmoltz instability or the Sedov- Taylor point-like explo- 
sion. Furthermore, the adoption of a variable softening 
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Fig. 49. Two components test, angular momenta conser- 
vations. Top panels: test TC3, bottom panels: test TC4. 
Shown are the total (solid lines), gas (dashed lines) and 
Dark Matter (dash-dotted lines) angular momenta as a 
function of time L{t)^ and conservations as a function of 
time L{t)/L(to) (same meaning of the symbols). Note that 
in the TC3 test the gas and Dark Matter lines are super- 
posed because they have identical angular momentum. 



lentghs algorithm allows for a higher degree of adaptivity 
without resulting in appreciable losses in terms of preci- 
sion and conservation of energy. While some typical flaws 
of SPH codes are still present (e.g., problems in the con- 
servation of vorticity because of spurious shear viscosity) , 
these new features clearly improve the method and pos- 
itively aid in curing well-known drawbacks of the SPH 
algorithm. 

It must be however pointed out that the new features 
must be extensively tested, especially in situations of in- 
terest for cosmological simulations (i.e., gas dynamics in 
presence of a Dark Matter and/or stellar components). 
Here, we have restricted our analysis to standard hydrody- 
namical problems. Other and new tests will be presented 
in the companion paper by Merlin et al. (2009, in prepa- 
ration), dedicated to the inclusion in EvoL of radiative 
cooling, chemical evolution, star formation, energy feed- 
back, and other non-standard algorithms. 



Acknowledgements. The authors would like to thank D.J. 
Price, S. Borgani, P.A. Thomas, L. Secco, S. Pasetto and G. 
Tormen for the useful discussions and for their help. All simula- 
tions were run using the Cineca facilities and/or the Monster 
cluster at the Department of Astronomy (Padova). The pic- 
tures have been produced using MatLab or Splash by D. 
Price IPricd (|200/t ). 



References 

Agertz, O., Moore, B., Stadel, J., Potter, D., Miniati, F., 
Read, J., Mayer, L., Gawryszczak, A., Kravtsov, A., 
Nordlund, A., Pearce, F., Quilis, V., Rudd, D., Springel, 
v.. Stone, J., Tasker, E., Teyssier, R., Wadsley, J., & 
Walder, R. 2007, MNRAS, 380, 963 

Alimi, J.-M., Serna, A., Pastor, C., & Bernabeu, G. 2003, 
Journal of Computational Physics, 192, 157 

Athanassoula, E., Fady, E., Lambert, J. C., & Bosma, A. 
2000, MNRAS, 314, 475 

Balsara, D. S. 1995, Journal of Computational Physics, 
121, 357 

Barnes, J. & Hut, P. 1986, Nat, 324, 446 

Bate, M. R. & Burkert, A. 1997, MNRAS, 288, 1060 

Benz, W. 1990, in Numerical Modelling of Nonlinear 

Stellar Pulsations Problems and Prospects, ed. J. R. 

Buchler, 269-+ 
Boss, A. P. & Bodenheimer, P. 1979, ApJ, 234, 289 
Brookshaw, L. 1986, Proceedings of the Astronomical 

Society of Australia, 6, 461 
Burkert, A. & Bodenheimer, P. 1993, MNRAS, 264, 798 
Garraro, G., Lia, C, & Chiosi, C. 1998, MNRAS, 297, 

1021 

Gouchman, H. M. P., Thomas, P. A., & Pearce, F. R. 1995, 

ApJ, 452, 797 
Dehnen, W. 2001, MNRAS, 324, 273 
Dolag, K., Vazza, F., Brunetti, G., & Tormen, G. 2005, 

MNRAS, 364, 753 
Dubinski, J. 1996, New Astronomy, 1, 133 
Einfeldt, B., Roe, P. L., Munz, C. D., & Sjogreen, B. 1991, 

Journal of Computational Physics, 92, 273 
Evrard, A. E. 1988, MNRAS, 235, 911 
Ewald, P. P. 1921, Annalen der Physik, 369, 253 
Fryxell, B., Olson, K., Ricker, P., Timmes, F. X., Zingale, 

M., Lamb, D. Q., MacNeice, P., Rosner, R., Truran, 

J. W., & Tufo, H. 2000, ApJS, 131, 273 
Gingold, R. A. & Monaghan, J. J. 1977, MNRAS, 181, 

375 

Gresho, P. M. & Chan, S. T. 1990, International Journal 

for Numerical Methods in Fluids, 11, 621 
Harfst, S., Theis, C, & Hensler, G. 2006, A&A, 449, 509 
Hernquist, L. 1993, ApJ, 404, 717 

Hernquist, L., Bouchet, F. R., & Suto, Y. 1991, ApJS, 75, 
231 

Hernquist, L. & Katz, N. 1989, ApJS, 70, 419 

Hu, X. Y. & Adams, N. A. 2005, Journal of Computational 

Physics, 213, 844 
Lia, C, Portinari, L., & Garraro, G. 2002, MNRAS, 330, 

821 

Liska, W. & Wendroff, B. 2003, SI AM Journal on Scientific 

Computing, 25, 995 
Lucy, L. B. 1977, AJ, 82, 1013 

McMillan, S. L. W. & Aarseth, S. J. 1993, ApJ, 414, 200 
Merlin, E. & Chiosi, C. 2007, A&A, 473, 733 
Merritt, D. 1996, AJ, 111, 2462 

Monaghan, J. J. 1988, Computer Physics 
Communications, 48, 89 



E. Merlin et al.: EvoL I 



33 



— . 1992, ARA&A, 30, 543 

— . 1997, Journal of Computational Physics, 136, 298 
— . 2000, Journal of Computational Physics, 159, 290 
— . 2002, MNRAS, 335, 843 
— . 2006, MNRAS, 365, 199 

Monaghan, J. J. & Gingold, R. A. 1983, Journal of 

Computational Physics, 52, 374 
Monaghan, J. J. & Lattanzio, J. C. 1985, A&A, 149, 135 
Morris, J. & Monaghan, J. J. 1997, Journal of 

Computational Physics, 136, 41 
Nelson, R. P. & Papaloizou, J. C. B. 1994, MNRAS, 270, 

1 

Noh, W. F. 1987, Journal of Computational Physics, 72, 
78 

Norman, M. L., Bryan, G. L., Harkness, R., Bordner, J., 
Reynolds, D., O'Shea, B., & Wagner, R. 2007, ArXiv 
e-prints 

Ott, F. & Schnetter, E. 2003, ArXiv Astrophysics e-prints, 

astro-ph/0303112 
Price, D. J. 2007, Publications of the Astronomical Society 

of Australia, 24, 159 
— . 2008, Journal of Computational Physics, 227, 10040 
Price, D. J. & Monaghan, J. J. 2005, MNRAS, 364, 384 
— . 2007, MNRAS, 374, 1347 
Romeo, A. B. 1998, A&A, 335, 922 
Rosswog, S. 2009, ArXiv e-prints 
Rosswog, S. & Price, D. 2007, MNRAS, 379, 915 
Saitoh, T. R. & Makino, J. 2008, ArXiv e-prints 
Salmon, J. K. & Warren, M. S. 1994, Journal of 

Computational Physics, 111, 136 
Schuessler, 1. & Schmitt, D. 1981, A&A, 97, 373 
Serna], A., Alimi, J.-M., & Chieze, J.-P. 1996, ApJ, 461, 

884 

Shirokov, A. 2007, ArXiv e-prints 

Sod, G. A. 1978, Journal of Computational Physics, 27, 1 
Springel, V. 2005, MNRAS, 364, 1105 
— . 2009, ArXiv Astrophysics e-prints, astro-ph/9506070 
Springel, V. & Hernquist, L. 2002, MNRAS, 333, 649 
Springel, V., Yoshida, N., & White, S. D. M. 2001, New 

Astronomy, 6, 79 
Steinmetz, M. & Mueller, E. 1993, A&A, 268, 391 
Thacker, R. J. & Couchman, H. M. P. 2000, ApJ, 545, 728 
Thomas, P. A. & Couchman, H. M. P. 1992, MNRAS, 257, 

11 

Wadsley, J. W., Stadel, J., & Quinn, T. 2004, New 

Astronomy, 9, 137 
Watkins, S. J., Bhattal, A. S., Francis, N., Turner, J. A., 

& Whitworth, A. P. 1996, A&AS, 119, 177 
Wetzstein, M., Nelson, A. F., Naab, T., & Burkert, A. 

2008, ArXiv e-prints 
Williams, P. R., Churches, D. K., & Nelson, A. H. 2004, 

ApJ, 607, 1 

Woodward, P. & Colella, P. 1984, Journal of 
Computational Physics, 54, 115 



Appendix A: A^-dimensional kernel 

The general form for the A/'-dimensional spline kernel func- 
tion (1 < AT < 3) is 



WN{r,e) =r]N X < 



(l-u'^^^u^ ifO<ix< 1, 
1(2 -u)^ ifl<ii<2, (A.l) 







if > 2, 



where 



r]N 



1 ifN = l, 

X < 15/(77r) if N = 2, 
^3/(27r) ifN = 3. 



(A.2) 



Derivatives are straightforwardly obtained from this ex- 
pression. The 1-D and 2-D kernels have been used where 
necessary in some of the hydrodynamical tests presented 
in Sect. [5l 

Appendix B: Summary of equations of motion 
and conservation 

Gravitational acceleration: 



dvi 



z,grav 



dt 



(B.l) 



i 



ri - Vi 



Ti dvi Tj dvj 



Hydrodynamical acceleration (for SPH particles only): 



dv 



i,hyd 



dt 



(B.2) 



Ci /m. 



ViWij{hj)+lIijViWij 



Specific internal energy evolution (for SPH particles only): 



dui 
'dt 



(B.3) 



E 



In these expressions. 



34 



E. Merlin et al.: EvoL I 



n* = 1 



dhi^dWijihi) 



dn, 



dhi 



(B.4) 



T,; 



drii ^ dei 



(B.5) 



^ _ ^ ^ d(j)ij{ei) 



drii 



da 



(B.6) 



and 



dWijihi) 



dhj ^ 

j 



(B.7) 



The smoothing kernel W is given by Eq. [H whereas the 
softened gravitational potential (j) is given by Eq. [HI 
The standard equation of state is that of an ideal gas, 



P = (7 - l)pu, 



(Bi 



where 7 is the adiabatic index (5/3 for a monatomic gas); a 
more general equation of state will be presented in Merlin 
et al. (2009, in preparation). 



